Documentation

General

Chicken, an important worldwide farmed animal, having the largest number and diverse morphology, not only serve as a good source of protein in our daily life but also used as a model organism for developmental biology, evolutionary genomics and basic medical research. Here, we developed a comprehensive chicken muti-omics database, Galbase, for providing five main functionalities: variation search, gene expression search, epigenomics search, phenotypic search, and a series of useful tools. In current version, Galbase contains 21,672,487 SNPs and 2,708,244 indels derived from 928 resequencing data of 47 domestic breeds and five subspecies of red jungle fowl. 429 transcriptome data and 379 epigenomic data were also integrated to identify gene expression pattern and cis-regulatory elements. We also collected 15,275 QTL entries and 7,526 associations manually curated from AnimalQTLdb, GWAS Atlas, and literature resources. Collectively, all these characteristics will make Galbase be a useful archive for in-depth analysis in biological research and breeding.

Related article:

Wang MS, Thakur M, Peng MS, Jiang Y, Frantz LAF, Li M, ... Han JL, Wu DD, Zhang YP. 863 genomes reveal the origin and domestication of chicken[J]. Cell research, 2020, 30(8): 693-701.

Ye S, Gao N, Zheng R, Chen Z, Teng J, Yuan X, Zhang H, Chen Z, Zhang X, Li J, Zhang Z. Strategies for obtaining and pruning imputed whole-genome sequence data for genomic prediction[J]. Frontiers in genetics, 2019, 10: 673.

Method

I. Pipeline of SNPs and indels calling

1. All raw sequencing data were obtained from Sequence Read Archive (SRA) of NCBI. We integrated a total of 928 chicken genome samples of 778 domestic chickens (including 47 breeds) and 150 Red jungle fowl (including five subspecies).
2. The adapter sequences and low-quality reads were trimmed using Trimmomatic v0.36. Then cleaned reads were mapped to the chicken reference genome GRCg6a (GCF_000002315.6) using the Burrows-Wheeler Aligner “BWA-MEM” with default parameters. The alignment BAM files were processed to sort reads, merge sample, and mark duplicates by Picard v2.1. Then, Genome Analysis Toolkit (GATK) HaplotypeCaller (version 3.7) were implemented to detect the SNPs in all samples. SNPs with only two alleles, missing rates < 0.1, and quality score ≥40 were extracted.
3. For indel calling, we first sifted structural variations for each sample by GATK with the SelectVariants-based method. Then, we applied the hard filter command "VariantFiltration" to exclude potential false-positive variant calls with the parameter –filterExpression "QD < 2.0 || FS > 200.0 || ReadPosRankSum < −20.0 || InbreedingCoeff < −0.8". Finally, we only retained the 1–30 bp indels for downstream analysis.
4. Annotation of SNPs and indels was performed by using snpEff.
5. Minor allele frequencies (MAF) for all chickens, and allele frequencies for each chicken group were calculated with PLINK and VCFtools.

II. Gene Expression

1. A total of 429 transcriptome data covering 44 tissues were downloaded from NCBI SRA.
2. The raw data were pre-processed to remove adaptor sequences and filter low-quality reads using Trimmomatic v0.36. Then high-quality reads were aligned against chicken GRCg6a genome with STAR (Version 2.5.1), and unmapped reads were further aligned by HISAT2 (Version 2.0.3-beta) to improve reads utilization.
3. We assembled transcripts and computed Transcripts Per Million (TPM) values for each genes by StringTie. And a tissue-specific index, Tau, was calculated by an in-house python script.

III. Epigenome Analysis

1. We downloaded 379 chicken epigenomics data, including histone ChIP-seq and ATAC-seq, from NCBI SRA.
2. Low-quality bases and adaptor sequences were trimmed by Trim Galore. Then, the cleaned reads were aligned against GRCg6a using Bowtie2 (v 2.2.8).
3. Low-quality and multiple-mapping reads were removed using Samtools with the option "-q 20". And Reads mapping to mitochondrial DNA were also removed.
4. Peak calling was performed by MACS (v.2.1.1), retaining only peaks with q-value < 0.05.

IV. Phenotype information

1. Phenotypic data were collected from Animal QTLdb, GWAS Atlas, and previously reported studies.
2. For Animal QTLdb, we retained entries with both chromosomal locations and traits in them.
3. Associations were collected from GWAS Atlas and literature resources.
4. The coordinates of all phenotypic data were converted to the GRCg6a genome version by using LiftOver.

V. liftOver chain file

1. We aligned the genome of GRCg7b to GRCg6a by minimap2, producing pairwise alignment result.
2. Then two utilities, maf-convert and axtChain, were used to produce liftOver chain files to facilitate users to coordinate conversion.

VI. Database implementation

1. High-quality SNPs, indels, gene expression, cis-regulatory elements, GWAS, and QTL, were processed with Perl scripts and stored in the MySQL database.
2. We use PHP Server Pages, HTML5, and JavaScript to implement search, data visualization and download.

Tutorial

Based on the standardized data processing and curation process, Galbase houses a wealth of omics data including reference genomes, annotations, high-quality genetic variants, transcriptomes, histone modifications, open chromatin regions, GWAS, and QTL. We divided these data into five main retrieval modules: (i) variation search, (ii) expression search, (iii) epigenomic search, (iv) phenotypic search, and (v) a series of useful tools. Each module has its own page, and features are linked to each other by gene symbol and chromosome location.

 

Fig 1. Data processing pipeline used to construct Galbase.

 

I. Variation search

The “Variation Module” has been designed to dynamically retrieve relevant information of SNPs and indels in tabular format, or in Gbrowse interface. Simply by specifying a variant rsID, a gene symbol, or a chromosomal location, users can easily obtain the results of a query for all annotated variation information, including chromosome, position, reference/alternative alleles, minor allele frequency (MAF), consequence type, variant ID, and allele frequency. Gbrowse is linked to this page with other sequence features. More filtering parameters can be set to obtain those SNPs that fall in different gene bodies, or non-protein-coding fraction. Moreover, a geographical map showing the allele frequency of five RJF subspecies and 47 chicken breeds can also be displayed, accompanied by each search. This function can help user to identify the breed or trait-associated variants.

  SNPs or indels search

Fig 2. Screenshots of a SNP data search and the results for examples.

II. Expression search

The “Expression Module” displays gene expression profiles of 429 RNA-seq data across 44 tissues in different stages, representing all the major organ systems of the chicken. In details, by searching for a specific gene symbol, Galbase provides three types of display formats of the expression data. The first one displays the expression matrix by heatmap, and also incorporate a tau value, which enables users easily find gene expression patterns and tissues with high expression abundance. The other two displays RNA-seq data in Gbrowse, such as the “expreBar” track, which generates a box-whisker plot showing the range of expression levels across samples; and the “RNAseqReadsCoverage” track, which displays normalized read coverage depth by converting BAM files to 1x sequencing depth. All these functions provide download links to expression matrix in CSV format or plain-text files.

 

Fig 3. Screenshots of chicken transcriptome information and the results for one example.

III. Epigenomic search

In “Epigenomic Module”, we collected four histone modification marks (H3K4me3, H3K27ac, H3K4me1, and H3K27me3) and one transcription factor, CCCTC-binding factor (CTCF), based on ChIP-seq, as well as one open chromatin marker, based on ATAC-seq to identify cis-regulatory elements. The epigenomics data can be accessed by table browser; both regulatory peaks and reads coverage (normalized 1x sequencing depth) can be visualized in Gbrowse; and the metadata, including sample information, experiment type, and epigenetic mark, can be displayed by heatmap views or typical ‘wiggle’ views in WashU Epigenome Browser. Tracks can be sorted, organized, dragged, and printed to PDF format based on user preferences.

 

Fig 4. Screenshots of chicken epigenomics data search and the results for one example.

IV. Phenotypic search

We collected 15,275 QTL entries and 7,526 associations manually curated from AnimalQTLdb, GWAS Atlas, and literature resources. Users can retrieve phenotypic data by different ways: search by gene symbol; search by variant rsid; find QTLs/associations by genome location; and find associated genes by a trait name or a keyword. We provide position search of two genome versions, including GRCg6a and GRCg7b genome assemblies.

 

Fig 5. Screenshots of chicken phenotypic data search.

V. A series of useful tools

1. Local UCSC genome browser

Users can search with a gene symbol, or a transcript name, or a genomic region to view SNPs, indels, gene expression, epigenomics signal, and phenotypic data in the global view. Currently, 1196 tracks for the GRCg6a assembly and 388 tracks for GRCg7b assembly have been released. The "PDF/PS" item under the "View" menu of navigation bar was used to generate a high quality image in PostScript or PDF formats.

2. Local WashU Epigenome Browser

Users can search with a gene symbol, or a transcript name, or a genomic region to view epigenomics signal pattern in our local WashU Epigenome Browser.

3. Alignment search tools (BLAT/BLAST)

We introduced two sequence alignment tools, webBlat and viroBLAST. The webBlat can be used to quickly search for homologous regions of a DNA or mRNA sequence, which can then be displayed in the browser. ViroBLAST can find regions of local similarity between sequences, which can be used to infer functional and evolutionary relationships between sequences.

4. Genome coordinate conversion tool (liftOver)

We also introduced a genome coordinate conversion tool, liftOver. The liftOver tool is used to translate genomic coordinates from one assembly version into another. Our database provides an online lift from GRCg6a to GRCg7b.

Project organizers