<jupyter_start><jupyter_text>Comparison between cattle graphs and human graphsWe conducted additional analysis to compare our cattle graphs with human graphs. Using sequence variants from 1000 human genomes project, we characterized graphs augmented with variants from four representative (and genetically distinct) human populations (YRI- African, GBR- European, STU- South Asia, and JPT- East Asia).<jupyter_code>#loading required library
library(tidyverse)
library(magrittr)
## We download data from github with this base url 
basepath <- "https://raw.githubusercontent.com/danangcrysnanto/bovine-graphs-mapping/master/part1_varselect/result/result_human"<jupyter_output><empty_output><jupyter_text>Plot relationship among samples with PCA<jupyter_code>## Looking at the four selected populations
## PCA performed with plink v1.9 ~2500 axis
## Plot the first and second axis of eigen vectors
datall   <- read.table(url(file.path(basepath,"all_1kg_pca.eigenvec")),header=TRUE,stringsAsFactors = FALSE)
head(datall)
## Calculate proportion explained from each axis from eigen value
### ratio should be based on the eigen val
egval  <- read.table(url(file.path(basepath,"all_1kg_pca.eigenval")),header=FALSE)
colnames(egval)  <- "valabs"
egval$vaprop  <- egval$valabs*100/sum(egval$valabs)
head(egval)<jupyter_output><empty_output><jupyter_text>First axis explained 1.77 % of the variations and second axis explained 0.56% from the total variation<jupyter_code>## select four population considered in this study
options(repr.plot.width=14, repr.plot.height=8)


datsel  <- datall  %>% filter(population %in% c("GBR","ITU","YRI","JPT"))
datsel$population  <- factor(datsel$population,levels=c("JPT","GBR","STU","YRI"),
                                       labels=c("JPT\nEast Asia","GBR\nEuropean","STU\nSouth Asia","YRI\nAfrican"))

pl1  <- ggplot()+
geom_point(data=datall,aes(x=pc1,y=pc2),col="grey80",size=1.5)+
geom_point(data=datsel,aes(x=pc1,y=pc2,col=population),size=2.5)+
theme_bw()+
theme(panel.border = element_rect(size=2),
      text=element_text(size=24),
      axis.title=element_text(face="bold"),
      panel.grid=element_blank(),
      legend.pos=c(0.5,0.1),
      legend.direction="horizontal")+
scale_colour_manual(values=c("violetred3","#0072B2", "#D55E00", "mediumpurple3"))+
labs(x="Principal component 1 (1.77%)",y="Principal component 2 (0.56%)",col="Population")

pl1<jupyter_output>Warning message:
“Removed 102 rows containing missing values (geom_point).”<jupyter_text>Plot nucleotide diversityNucleotide diversity (phi) calculated based on the vcftools with 10 kb non-overlapped windows<jupyter_code>datdiv <- read.table(url(file.path(basepath,"phi_human.tsv")),header=FALSE,stringsAsFactors = FALSE)
colnames(datdiv) <- c("phi","pop")
head(datdiv)
datsum  <- datdiv %>% group_by(pop)  %>% summarise(nc=mean(phi))
datdiv$pop <- factor(datdiv$pop,levels=c("JPT","GBR","STU","YRI"))
pl2  <- ggplot(datdiv,aes(x=pop,y=phi))+
  geom_boxplot(aes(col=pop),show.legend=FALSE)+
  geom_text(data=datsum,aes(x=pop,y=-0.0005,label=round(nc,5)), inherit.aes = FALSE,size=6)+
  theme_bw()+
  theme(panel.border = element_rect(size=2),
        text=element_text(size=24),
        axis.title=element_text(face="bold"),
        panel.grid = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(size=20))+
  scale_y_continuous(limits=c(-0.0005,0.015))+
  scale_color_manual(values=c("violetred3","#0072B2", "#D55E00", "mediumpurple3"))+
  labs(x="Populations",y="Nucleotide diversity")

pl2<jupyter_output><empty_output><jupyter_text>Plot graph mapping performance We constructed series of graphs from four human population with filter on allele frequency. We then simulated 10 millions reads from selected samples from each population. Then, we map simulated reads to all graphs from corresponding population. We calculated and plotted the mapping error as indicated below.<jupyter_code># the data consists of 4 columns with population, af threshold, total incorrect mapping, and mapping type
## subset to include only paired end
dathum <- read.table(url(file.path(basepath,"dathum_freq.tsv")),stringsAsFactors = FALSE,header=TRUE)
dathum_pe  <- dathum  %>% select(humpop,afcor,pecor2)
dathum_pe$maptype  <- "pe"
colnames(dathum_pe)  <- c("humpop","af","incor2","maptype")
head(dathum_pe)
dathum_pe$humpop <- factor(dathum_pe$humpop,levels=c("JPT","GBR","STU","YRI"))

pl3  <- ggplot(dathum_pe,aes(x=af,y=incor2,col=humpop))+
  geom_point(size=3)+
  geom_line()+
  theme_bw()+
  scale_x_continuous(breaks = c(0,0.05,seq(0.1,1.1,0.1))) +
  theme(panel.border = element_rect(size=2),
        text=element_text(size=24),
        axis.title=element_text(face="bold"),
        axis.text.x=element_text(size=14),
        panel.grid=element_blank(),
        strip.background = element_blank(),
        strip.text = element_blank(),
        legend.position = "bottom")+
  labs(x="Allele frequency threshold",y="Mapping error (%)",col="Population")+
  scale_color_manual(values=c("violetred3","#0072B2", "#D55E00", "mediumpurple3"))

pl3<jupyter_output><empty_output><jupyter_text>Plot of the TPR-FPR from the best graph from each populationWe selected the best graph from each population and calculate TPR and FPR based on the mapping quality<jupyter_code># The table divided by allele frequency bin and corresponding TPR-FPR from each bin
datroc <- read.table(url(file.path(basepath,"best_roc_human.tsv")),header=TRUE,stringsAsFactors = FALSE)
datroc_pe  <- datroc  %>% filter(mode=="pe")
head(datroc_pe)
datroc_pe$population <- factor(datroc_pe$population,levels=c("JPT","GBR","STU","YRI"))

pl4 <- ggplot(datroc_pe,aes(x=FPR/100,y=TPR/100,col=population))+
  geom_point(size=3)+
  geom_line()+
  scale_x_log10(breaks=c(10^-6,10^-5,10^-4,10^-3,10^-2,10^-1,0),
                labels=scales::trans_format("log10",scales::math_format(10^.x)))+
  theme_bw()+
  scale_y_continuous(breaks=c(0.984,0.988,0.992,0.996,1.000),limits=c(0.984,1.000))+
  theme(panel.border = element_rect(size=2),
        text=element_text(size=24),
        axis.title=element_text(face="bold"),
        panel.grid=element_blank(),
        strip.background = element_blank(),
        strip.text = element_blank(),
        legend.position = "bottom")+
  scale_color_manual(values=c("violetred3","#0072B2", "#D55E00", "mediumpurple3"))+
  guides(col = guide_legend(nrow = 1))+
  labs(x="FPR",y="TPR",col="Population")

pl4<jupyter_output><empty_output><jupyter_text>Comparison with random variant selectionWe selected the random variants as the number of variants filtered based on frequency and look into the mapping profile.<jupyter_code>##The data consists of breed, allele frequency threshold, incorrect mapping, and number of variants at given bin
datrand <- read.table(url(file.path(basepath,"random_af_human.tsv")),header = TRUE, stringsAsFactors = FALSE)
datsum <- datrand %>% group_by(af,maptype) %>% summarise(nc_pe=mean(pecor),
                                                 nc_se=mean(secor),
                                                 totvar=mean(totvar))
head(datrand)
head(datsum)
pl5  <-  ggplot(datrand,aes(x=totvar,y=pecor*10^2/10^7))+
  geom_point(aes(col=population,shape=maptype),size=2,alpha=0.7)+
  geom_line(data=datsum,aes(y=nc_pe*10^2/10^7,linetype=maptype),size=1)+
  theme_bw()+
  scale_color_manual(values=c("violetred3","#0072B2", "#D55E00", "mediumpurple3"))+
  theme(panel.border = element_rect(size=2),
        text = element_text(size=24),
        panel.grid = element_blank(),
        axis.title = element_text(face="bold"),
       legend.position="bottom")+
  scale_x_continuous(labels = c(0,100,200,300), limits = c(0,3*10^5))+
  scale_shape_manual(values=c(17,16)) +
  labs(x="Number of variants (x1000)",y="Mapping error (%)")

pl5<jupyter_output><empty_output><jupyter_text>Plotting the single-end results<jupyter_code>##af plot as in figure (pl3)
dathum_se  <- dathum  %>% select(humpop,afcor,secor2)
dathum_se$maptype  <- "se"
colnames(dathum_se)  <- c("humpop","af","incor2","maptype")
dathum_se$humpop <- factor(dathum_se$humpop,levels=c("JPT","GBR","STU","YRI"))

pl3_se  <- ggplot(dathum_se,aes(x=af,y=incor2,col=humpop))+
  geom_point(size=3)+
  geom_line()+
  theme_bw()+
  scale_x_continuous(breaks = c(0,0.05,seq(0.1,1.1,0.1))) +
  theme(panel.border = element_rect(size=2),
        text=element_text(size=24),
        axis.title=element_text(face="bold"),
        axis.text.x=element_text(size=14),
        panel.grid=element_blank(),
        strip.background = element_blank(),
        strip.text = element_blank(),
        legend.position = "bottom")+
  labs(x="Allele frequency threshold",y="Mapping error (%)")+
  scale_color_manual(values=c("violetred3","#0072B2", "#D55E00", "mediumpurple3"))

pl3_se
#tpr-fpr as in figure (pl) 4
datroc_se  <- datroc  %>% filter(mode=="se")
datroc_se$population <- factor(datroc_pe$population,levels=c("JPT","GBR","STU","YRI"))

pl4 <- ggplot(datroc_se,aes(x=FPR/100,y=TPR/100,col=population))+
  geom_point(size=3)+
  geom_line()+
  scale_x_log10(breaks=c(10^-6,10^-5,10^-4,10^-3,10^-2,10^-1,0),
                labels=scales::trans_format("log10",scales::math_format(10^.x)))+
  theme_bw()+
  theme(panel.border = element_rect(size=2),
        text=element_text(size=24),
        axis.title=element_text(face="bold"),
        panel.grid=element_blank(),
        strip.background = element_blank(),
        strip.text = element_blank(),
        legend.position = "bottom")+
  scale_color_manual(values=c("violetred3","#0072B2", "#D55E00", "mediumpurple3"))+
  guides(col = guide_legend(nrow = 1))+
  labs(x="FPR",y="TPR",col="Population")

pl4
### random sampling plot as in figure (pl) 5
datrand <- read.table(url(file.path(basepath,"random_af_human.tsv")),header = TRUE, stringsAsFactors = FALSE)
datsum <- datrand %>% group_by(af,maptype) %>% summarise(nc_pe=mean(pecor),
                                                 nc_se=mean(secor),
                                                 totvar=mean(totvar))
pl5_se  <-  ggplot(datrand,aes(x=totvar,y=secor*10^2/10^7))+
  geom_point(aes(col=population,shape=maptype),size=2,alpha=0.7)+
  geom_line(data=datsum,aes(y=nc_se*10^2/10^7,linetype=maptype),size=1)+
  theme_bw()+
  scale_color_manual(values=c("violetred3","#0072B2", "#D55E00", "mediumpurple3"))+
  theme(panel.border = element_rect(size=2),
        text = element_text(size=24),
        panel.grid = element_blank(),
        axis.title = element_text(face="bold"),
        legend.position="bottom")+
  scale_x_continuous(labels = c(0,100,200,300), limits = c(0,3*10^5))+
  scale_shape_manual(values=c(17,16)) +
  labs(x="Number of variants (x1000)",y="Mapping error (%)")

pl5_se
sessionInfo()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15892193308550187, 'fraction_numerical': 0.0379182156133829, 'mean_word_length': 5.318849089841456, 'pattern_counts': {'":': 0, '<': 69, '<?xml version=': 0, '>': 43, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '620430', 'license': 'MIT', 'n_tokens_mistral': 4168, 'n_tokens_neox': 3861, 'n_words': 633, 'path': 'part1_varselect/analysis/part1_human_bovine_comparison.ipynb', 'repo_name': 'danangcrysnanto/bovine-graphs-mapping', 'size': '10760'}