#An example of how to import a QTL annotation file#qtl.inp <- import_gff_gtf(db_file="QTL_db.gff",file_type="gff")#An example of how to import a gene annotation file#qtf.inp <- import_gff_gtf(db_file="Gene_db.gtf",file_type="gtf")

Gene annotation

For example, let`s run a gene and QTL annotation using the QTLwindowsdataset without additional intervals (upstream and downstream, using theinterval=0 argument) from the windows coordinates:

#Running gene annotationout.genes<-find_genes_qtls_around_markers(db_file=gtfGenes,marker_file=QTLmarkers,method ="gene",marker ="snp",interval =500000,nThreads =1)

You are using the method: gene with snp

## Warning: executing %dopar% sequentially: no parallel backend registered

#Checking the first rows from the output fileDT::datatable(out.genes,rownames =FALSE,extensions ='FixedColumns',options =list(scrollX =TRUE))

#Checking the dimensions of the output filedim(out.genes)

[1] 652 15

The gene annotation resulted in 652 genes within the 1 Mb interval(500 Kb upstream and 500 Kb downstream) from the candidate markers.

QTL annotation

#Running QTL annotationout.qtls<-find_genes_qtls_around_markers(db_file=gffQTLs,marker_file=QTLmarkers,method ="qtl",marker ="snp",interval =500000,nThreads =1)

You are using the method: qtl with snp Starting QTL searching using5e+05 bp as interval Preparing output file for QTL annotation

#Checking the first rows from the output fileDT::datatable(out.qtls,rownames =FALSE,extensions ='FixedColumns',options =list(scrollX =TRUE))

#Checking the dimensions of the output filedim(out.qtls)

[1] 3020 37

The QTL annotation resulted in 3020 QTLs within the 1 Mb interval(500 Kb upstream and 500 Kb downstream) from the candidate markers.

HINT:It is important to highlight that bothoutputs are composed by all the columns in the input file (marker_file)plus the annotation columns from the gtf or gff file used in theannotation processes. Additionally, each row from the input file isrepeated as many times as an annotation (gene or QTL) record wasidentified in the determined interval. Consequently, for example, for aninput file composed of three genomic coordinates where 4 genes areannotated in the interval determined by the user, the output file offind_genes_qtls_around_markers() will contain 12 rows.

This output can be easily handled by summary functions in R, such astable(), to obtain information such as the total number of genes andQTLs, the number of genes and QTLs annotated per variants, etc. The sameoutput file generated can be used as an input file for the other set ofGALLO functions. This additional set of functions allows the graphicalvisualization of the results obtained including the summary of QTL typesand traits annotated, the overlapping among traits, populations andother kinds of groups, as well as the relationship between candidategenetic regions and the annotated genes and QTLs.

Checking overlap among groups

• Arguments from overlapping_among_groups

overlapping_among_groups(file, x, y)

file: A dataframe with the data and groupingfactor

x: The grouping factor to be compared

y: The data to be compared among the levels of thegrouping factor

#Removing the duplicated records for Gene IDout.genes.unique<-out.genes[!duplicated(out.genes[c("Reference","gene_id")]),]#Calculating the number of shared genes between different studiesoverlap.genes<-overlapping_among_groups(file=out.genes.unique,x="Reference",y="gene_id")overlap.genes

## $N##                         Feugang et al. (2010) Buzanskas et al. (2017)## Feugang et al. (2010)                     473                       1## Buzanskas et al. (2017)                     1                      42## ## $percentage##                         Feugang et al. (2010) Buzanskas et al. (2017)## Feugang et al. (2010)                    1.00                       0## Buzanskas et al. (2017)                  0.02                       1## ## $combined##                         Feugang et al. (2010) Buzanskas et al. (2017)## Feugang et al. (2010)   "473 (1)"             "1 (0)"                ## Buzanskas et al. (2017) "1 (0.02)"            "42 (1)"

The output of this function is composed by a list with threematrices: 1) A matrix with the number of overllaping data; 2) A matrixwith the percentage of overlapping; 3) A matrix with the combination ofthe two previous ones.

Now it is possible to plot the overlapping among the studies usingthe plot_overlapping function.

• plot_overlapping

plot_overlapping(overlapping_matrix, nmatrix, ntext, group,labelcex = 1)

overlapping_matrix: The list obtained inoverlapping_amoung_groups function

nmatrix: An interger from 1 to 3 indicating wichmatrix will be used to plot the overlapping, where: 1) A matrix with thenumber of shared data; 2) A matrix with the percentage of overlapping;3) A matrix with the combination of the two previous matrices

ntext: An interger from 1 to 3 indicating wichmatrix will be used as the text matrix for the heatmap, where: 1) Amatrix with the number of shared data; 2) A matrix with the percentageof overlapping; 3) A matrix with the combination of the two previousone

group: A vector with the size of groups. This vectorwill be plotted as row and column names in the heatmap

labelcex: A numeric value indicating the size of therow and column labels

#Creating grouping labels (names of the references)group.labels<-unique(out.genes.unique$Reference)

#plotting the overlapping informationplot_overlapping(overlap.genes,nmatrix=2,ntext=3,group=group.labels,labelcex =1)

Figure 1: Gene overlapping among studies.

HINT:The output matrices are not necessarilysymmetric. For example, the percentage of shared genes between Feugantet al. (2010) and Buzanskas et al. (2017) is not the same percentage ofshared genes between Buzanskas et al. (2017) and Feugang et al. (2010).The percentage of shared genes is calculated based on the total numberof genes from the grouping factor “A” that are also present in the group“B” and vice-versa. Once the total number of genes in each group can bedifferent, the percentage of genes (or any other data) from A sharedwith B can not be the same as B shared with A. The diagonals in thematrices represent the total number of genes in each study.

Plot QTL annotation

The GALLO package also allows the user to plot the QTL informationannotated using the plot_qtl_info function. This function can create apie plot (using the option “qtl_type” in the argument qtl_plot) with thepercentage of QTL classes and/or bar using the option “qtl_name” in theargument qtl_plot) plots with the percentage of each trait that composedthe slice of the pie plot for a respective QTL class.

• plot_qtl_info

  plot_qtl_info(qtl_file, qtl_plot = c(“qtl_type”, “qtl_name”),n =“all”, qtl_class = NULL, …)

  qtl_file: The output fromfind_genes_qtls_around_markers function

  qtl_plot: “qtl_type” or”qtl_name”

  n: Number of QTLs to be plotted when the qtl_nameoption is selected

  qtl_class: Class of QTLs to be plotted when theqtl_name option is selected

#plotting the percentage of each QTL class annoattedoldpar<-par(mar=c(0.5,15,0.5,1))plot_qtl_info(out.qtls,qtl_plot ="qtl_type",cex=1.5)

par(oldpar)

The Reproduction QTL class was the most frequent annotated QTL class(44.01%), followed by Milk (31.06%), Production (10.1%), Exterior(8.21%), Meat and Carcass (4.97%) and health (2.55%). The percentage ofeach QTL classe can give us a better idea about how frequently ourcandidate regions are associated with each kind of QTL class. However,the investigation bias (number of studies) for a QTL traits, such asMilk related QTLs, will result in much more QTLs annotated for thattrait. Consequently, we need to investigate in deep the relationshipbetween our candidate regions and the previous reported QTLs. Thisanalysis can help us to identify specialized regions across the genomefor a trait and/or pleiotropic and espistatic effects in our candidateregions.

Using the same fucntion, plot_qtl_info, now we can check therespective percentages of each trait that was annotated and composedthat QTL class. In order to do this, we need to inform the“Reproduction” class for the qtl_class argument.

#Setting margin parameter to better fit the axis labelsoldpar<-par(mar=c(5,20,1,1))#plotting the percentage of each trait annoatted as a Reproduction QTLplot_qtl_info(out.qtls,qtl_plot ="qtl_name",qtl_class="Reproduction")

Figure 2: Percentage of each trait annotated as aReproduction QTL in the candidate intervals.

par(oldpar)

This plot allows to check the most representative traits within eachkind of QTL classes.

HINT:The user can easily include this functionin a for loop and plot the percentages for all the QTL classes.Additionally, the number of traits exhibited in the plots can be passedto the argument n. This can be performed using the following code, forexample:

#Plotting percentage of the top 10 most frequent traits in all QTL classes#(This is just an example code, the user do not need to execute#this command for this tutorial)QTL_classes<-unique(out.qtls$QTL_type)for(cin QTL_classes){  tmp.file.name<-paste(c,".png",sep="")png(tmp.file.name,w=1500,h=900)plot_qtl_info(out.qtls,qtl_plot ="qtl_name",qtl_class=c,n=10)dev.off()}

qtl_enrich

qtl_enrich(qtl_db, qtl_file, qtl_type = c(“QTL_type”,“trait”),enrich_type = c(“genome”, “chromosome”), chr.subset = NULL,nThreads = NULL, padj = c(“holm”, “hochberg”,“hommel”, “bonferroni”,“BH”, “BY”, “fdr”, “none”))

qtl_db: The .gff file that can be downloaded fromAnimal QTlLdb

qtl_file: The output fromfind_genes_qtls_around_markers function

qtl_type: A character indicating which type ofenrichment will be performed. “QTL_type” indicates that the enrichmentprocesses will be performed for the QTL classes, while “Name” indicatesthat the enrichment analysis will be performed for each traitindividually.

enrich_type: A character indicating if theenrichment analysis will be performed for all the chromosomes (“genome”)or for a subset of chromosomes (“chromosome”). If the “genome” option isselected, the results reported are the merge of all chromosomes.

chr.subset: If enrich_type is equal “chromosome”, itis possible to define a subset of chromosomes to be analyzed. Thedefault is equal NULL. Therefore, all the chromosomes will beanalyzed.

nThreads: The number of threads used.

padj: The alogorithm for multiple testing correctionto be adopted (“holm”, “hochberg”, “hommel”, “bonferroni”, “BH”,“BY”,“fdr”, “none”).

As an example, we are going to perform a enrichment analysis for allthe QTL information annotated around the candidate markers using achromosome-based enrichment analysis. The adjusted p-values will becalculated based on False-Discovery Rate (FDR).

This step might take some minutes to run depending of the user`ssystem.

#QTL enrichment analysisout.enrich<-qtl_enrich(qtl_db=gffQTLs,qtl_file=out.qtls,qtl_type ="Name",enrich_type ="chromosome",chr.subset =NULL,padj ="fdr",nThreads =1)

Staring QTL enrichment analysis for QTL class

## End of QTL enrichment analysis

#Checking the enriched QTLsDT::datatable(out.enrich[order(out.enrich$pvalue),],rownames =FALSE,extensions ='FixedColumns',options =list(scrollX =TRUE))

The output from the qtl_enrich function is a data frame composed by 9columns: 1) QTL: The QTL class or trait used for the enrichment; 2) CHR:The chromosome for that specific QTL or trait (if the option“chromosome” is informed to the argument enrich_type); 3) N_QTLs: Numberof observed QTLs or traits in the dataset; 4) N_QTLs_db: Number of eachannotated QTL in the qTL database; 5) Total_annotated_QTLs: Total numberof annotated QTLs; 6) Total_QTLs_db: Total number of QTLs in the QTLdatabase; 7) pvalue: P-value for the enrichment analysis; 8) adj.pval:The adjusted p-value based on the multiple test correction selected bythe user; 9) QTL_type= The QTL type for each annotated trait.

HINT:It is possible to match the traits in theqtl_enrich output with the respective QTL classes from the QTLannotation output. Consequently, the enrichment results can be filteredby the QTL classes, such as Reproduction, Meat and Carcass, Production,etc.

QTLenrich_plot

QTLenrich_plot(qtl_enrich, x, pval)

qtl_enrich: The output from qtl_enrich function

x: ID column to be used from the qtl_enrichoutput

pval: P-value to be used in the plot. If “p_value”informed, a non-adjusted pvalue will be plotted. If “p.adj” informed,the adjusted p-value from the qtl enrichment analysis will be plotted.

Before plot the enrichment results, a new ID column will be createdin order to make easier to identify the enrichment results perchromosome. Additionally, we are going to match the QTL class for eachtrait and filter the top 10 enriched QTLs.

#Creating a new ID composed by the trait and the chromosomeout.enrich$ID<-paste(out.enrich$QTL," - ","CHR",out.enrich$CHR,sep="")#Match the QTL classes and filtering the Reproduction related QTLsout.enrich.filtered<-out.enrich[which(out.enrich$adj.pval<0.05),]#Here we are going to truncate the FDR values at -log10(5^-50) in order to provide a better visualizationout.enrich.filtered$new_pval<-out.enrich.filtered$adj.pvalout.enrich.filtered[which(out.enrich.filtered$new_pval<(5^-50)),"new_pval"]<-(5^-50)

#Plotting the enrichment results for the QTL enrichment analysisQTLenrich_plot(out.enrich.filtered,x="ID",pval="new_pval")

Figure 3: Top 10 enriched traits identified in theQTL enrichment analysis. The area of the bubbles represents the numberof observed QTLs for that class, while the color represents the p-valuescale (the darker the color, smaller the p-value). Additionally, thex-axis shows the richness factor for each QTL, representing the ratio ofnumber of QTLs and the expected number of that QTL.

As any enrichment analysis, the quality of annotation, in thiscase, the number of association studies for a specific trait, candirectly affect the enrichment results.

relationship_plot

relationship_plot(qtl_file, x, y, grid.col = “gray60”, degree =90,canvas.xlim = c(-2, 2), canvas.ylim = c(-2, 2), cex)

qtl_file: The output fromfind_genes_qtls_around_markers function

x: The first grouping factor, to be plotted in theleft hand side of the chord plot

y: The second grouping factor, to be plotted in theleft hand side of the chord plot

grid.col: A character with the grid color for thechord plot or a vector with different colors to be used in the gridcolors. Note that when a color vector is provided, the lenght of thisvector must be equal the number of sectors in the chord plot

degree: A numeric value corresponding to thestarting degree from which the circle begins to draw. Note this degreeis always reverse-clockwise

canvas.xlim: The coordinate for the canvas in thex-axis. By default is c(-1,1)

canvas.ylim: The coordinate for the canvas in they-axis. By default is c(-1,1)

cex: The size of the labels to be printed in theplot

#Filtering the output from QTL annotation#Creating a new ID to filter the top 10 enriched QTLsout.qtls$ID<-paste(out.qtls$Name," - ","CHR",out.qtls$CHR,sep="")out.enrich.filtered<-out.enrich.filtered[order(out.enrich.filtered$adj.pval),]out.qtls.filtered<-out.qtls[which(out.qtls$ID%in%out.enrich.filtered$ID[1:10]),]#Creating color scheme based on the Referencesout.qtls.filtered[which(out.qtls.filtered$Reference=="Feugang et al. (2010)"),"color_ref"]<-"purple"out.qtls.filtered[which(out.qtls.filtered$Reference=="Buzanskas et al. (2017)"),"color_ref"]<-"pink"#Creating a color vector filled with black for all the traits abbreviation#and with the respective colors for each referencecolor.grid<-c(rep("black",length(unique(out.qtls.filtered$Abbrev))),out.qtls.filtered[!duplicated(out.qtls.filtered$SNP.reference),"color_ref"])#Naming the vectornames(color.grid)<-c(unique(out.qtls.filtered$Abbrev),unique(out.qtls.filtered$SNP.reference))#Plotting the relationship plot using the grid color created aboverelationship_plot(qtl_file=out.qtls.filtered,x="Abbrev",y="SNP.reference",cex=1,gap=2.5,degree =60,canvas.xlim =c(-1.5,1.5),canvas.ylim =c(-1.5,1.5),grid.col = color.grid)

Figure 4: Chord plot showing the relationship betweenthe positional candidate markers identified in the GWAS (right-handside) and the enriched QTLs (abbreviations in the left -hand side). Thecolors correspond to each study, where markers identified by Feugang etal. (2010) are shown in purple, while the markers identified byBuzanskas et al. (2017) are shown in pink.

The GALLO package provides a user-friendly and straightforwardenvironment to perform gene and QTL annotation, visualization, datacomparison and QTL enrichment for functional studies in livestockspecies. Consequently, the use of GALLO in the analysis of datagenerated from high-throughput methodologies may improve theidentification of hidden pattern across datasets, datamining of previousreported associations, as well as the efficiency in the scrutinizationof the genetic architecture of complex traits in livestock.

This was a short example of GALLO usage. The user can perform severaldifferent approaches in order to explore GALLO potential.