Documentation

General

The Ruminant Genome Database (RGD) is a comprehensive repository of integrated genomics, transcriptomics, and epigenetics data of Ruminantia, a suborder which includes economically important livestock: cattle, yak, sheep, and goats. The RGD builds a feature-rich UCSC Genome Browser to the local server, and implements open-source software and MySQL databases providing access to all the data contributed by our Ruminant Genome Project (RGP) and other public database. The current version of the RGD contains data of 55 ruminant genomes, whole genome alignments of 55 ruminants and 12 mammals, five scales of conservative evaluation, 1151 samples of RNA-seq data, as well as regulatory regions predicted by utilizing 32 ruminant epigenomic and 726 human epigenomic datasets. The gene ontology (GO) terms, KEGG pathway, and gene expression heatmap can be queried with gene symbols. RGD also provides online alignment tools such as webBlat and NCBI’s wwwBLAST server against the goat ARS1 assembly. All the data were integrated into UCSC Genome Browser to facilitate global presentation. Collectively, all these features make the RGD a useful archive for in-depth analysis in biology and evolution of ruminants.

Methods

1. Cross-species annotation

Firstly, we aligned goat whole genome (GCA_001704415.1) to the sheep (GCA_000298735.2), cattle (GCA_000003055.3 and GCA_000003055.5), and human (GCA_000001405.22) genome respectively to generate pairwise alignments using LAST v88555. Then gene coordinates of the sheep, cattle and human were converted to goat genome using LiftOver tool to improve goat gene annotation.

step 1: lastdb -uNEAR -cR11 ref.fa.db ref.fa
step 2: lastal -P20 -m100 -E0.05 ref.fa.db query.fa | last-split -m1 > query.maf
step 3: maf-swap query.maf | last-split | maf-swap | last-split | maf-sort > query.LAST.maf
step 4: maf-convert psl query.LAST.maf > query.psl
step 5: axtChain -linearGap=loose -psl query.psl ref.2bit query.2bit queryToRef.chain
step 6: liftOver -minMatch=0.9 -genePred ref.gene.genePred queryToRef.chain query.gene.map query.gene.unMapped (for sheep and cattle)
        liftOver -minMatch=0.2 -minBlocks=0.5 -genePred ref.gene.genePred queryToRef.chain query.gene.map query.gene.unMapped (for human)

2. 67-way multiz align

We built whole-genome alignments for 67 species (55 ruminants and 12 mammals) with LAST v88555 and MULTIZ, using the goat ARS1 assembly as the reference.

step 1: lastdb -uNEAR -cR11 ref.fa.db ref.fa
step 2: lastal -P20 -m100 -E0.05 ref.fa.db query.fa | last-split -m1 > query.maf
step 3: maf-swap query.maf | last-split | maf-swap | last-split | maf-sort > query.LAST.maf
step 4: perl maf.rename.species.S.pl query.LAST.maf refSpeciesName querySpeciesName query.Final.maf > query.stat
step 5: multiz query1.Final.maf query2.Final.maf 0 all > query1.query2.mutiz.maf

Get perl script: maf.rename.species.S.pl

3. Conservation scores

Both phastCons and phyloP conservation scores have been computed separately for three groups of organisms: Bovidae (38 species), Pecora (53 species) and Ruminatia (54 species).

[generate non-conservative and conservative models]

step 1: phyloFit --tree tree.nwk --msa-format FASTA --out-root nonconserved-4d.mod 4DTV.fasta
step 2: msa_split chr.maf --in-format MAF --windows 100000,0 --out-root split_maf/chr --out-format SS --min-informative 1000 --between-blocks 5000
step 3: phastCons --estimate-rho window.ss --no-post-probs window.ss nonconserved-4d.mod.mod
step 4: phyloBoot --read-mods window.cons.list --output-average Average.chr.cons.mod
        window.cons.list: 6.10000001-10100000.ss.cons.mod,6.1000001-1100000.ss.cons.mod,6.100001-200000.ss.cons.mod,6.100099918-100200000.ss.cons.mod,6.100200001-100300000.ss.cons.mod ...

[ phastCons ]

phastCons --most-conserved phastCons.chr.bed --score chr.maf Average.chr.cons.mod,nonconserved-4d.mod.mod > phastCons.chr.wig

[ phyloP ]

phyloP --wig-scores --mode CONACC --method LRT nonconserved-4d.mod.mod chr.maf > phyloP.chr.wig

4. RNA-seq analysis

Pipeline: RNA-seq data with a reference genome

step 1: [Quality_control]
        Single-read: java -Xmx30g -jar trimmomatic-0.36.jar SE -phred33 -threads 10 sample.fq.gz sample.clean.fq.gz LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:40 TOPHRED33 > sample.log
        Paired-end: java -Xmx30g -jar trimmomatic-0.36.jar PE -phred33 -threads 10 sample_1.fq.gz sample_2.fq.gz sample_1.clean.fq.gz sample_1.unpaired.fq.gz sample_2.clean.fq.gz sample_2.unpaired.fq.gz LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:40 TOPHRED33 > sample.log
step 2: [STAR_index]
        STAR --runMode genomeGenerate --genomeDir index --genomeFastaFiles ref.fa --runThreadN 24 --sjdbGTFfile ref.gtf --sjdbOverhang 149
step 3: [Hisat2_index]
        hisat2-build ref.fa ref.hisat2.index
step 4: [Splice_sites]
        python3.5 hisat2_extract_splice_sites.py ref.gtf > ref.hisat.splicing.site
step 5: [Mapping]
        STAR --genomeDir $INDEX --readFilesIn $IN/${i}_1.clean.fq.gz $IN/${i}_2.clean.fq.gz --readFilesCommand zcat --runThreadN 8 --outFilterMultimapNmax 1 --outFilterIntronMotifs RemoveNoncanonical Unannotated --outFilterMismatchNmax 10 --outSAMstrandField intronMotif --outSJfilterReads Unique --outSAMtype BAM Unsorted --outReadsUnmapped Fastx --outFileNamePrefix $OUT/${i}
        hisat2 --dta-cufflinks --no-mixed --no-discordant -p 10 --known-splicesite-infile ref.hisat.splicing.site -x ref.hisat2.index -1 $OUT/${i}Unmapped.out.mate1 -2 $OUT/${i}Unmapped.out.mate2 -S $OUT/${i}.sam
        java -Djava.io.tmpdir=$TMP -jar picard.jar CleanSam I=$OUT/${i}Aligned.out.bam O=$OUT/${i}.STAR.bam
        java -Djava.io.tmpdir=$TMP -jar picard.jar CleanSam I=$OUT/${i}.sam O=$OUT/${i}.hisat.bam
        java -Djava.io.tmpdir=$TMP -jar picard.jar MergeSamFiles I=$OUT/${i}.STAR.bam I=$OUT/${i}.hisat.bam SORT_ORDER=coordinate O=$BAM/${i}.bam
step 6: [Assemble_transcripts]
        stringtie sample.bam -G ref.gff -c 3 -j 2 -f 0.05 -o sample.gtf -p 16
        stringtie --merge -G ref.gff -o merge.gff allSample.gtf.list -g 100 -m 200 -c 0.8
        stringtie -e -B -p 8 -G merge.gff -o sample.gtf sample.bam
step 7: FPKMs (Fragments Per Kilobase per Million mapped reads) of these transcripts and genes in each sample were computed by Ballgown (Version 2.2.0)

5. ChIP-seq and ATAC-seq analysis

Pipeline: ChIP-seq data processing

step 1: [BWA mappping]
        Index: bwa index ref.fa
        Single-read: bwa mem ref.fa -t 8 -M -R '@RG\tID:sample\tLB:sample\tPL:ILLUMINA\tSM:sample' sample.fastq | samtools view -bS -@ 8 > sample.bam
        Paired-end: bwa mem ref.fa -t 8 -M -R '@RG\tID:sample\tLB:sample\tPL:ILLUMINA\tSM:sample' sample_1.fastq sample_2.fastq | samtools view -bS -@ 8 > sample.bam
step 2: [Merge, sort, and filter bam]
        java -Djava.io.tmpdir=/tmp/ -jar picard.jar MergeSamFiles I=sample_line1.bam I=sample_line2.bam I=sample_line3.bam SORT_ORDER=coordinate O=sample.bam
        samtools sort -@ 8 -o sample.sorted.bam -T sample.sorted.tmp -O bam sample.bam
        samtools view -q 20 sample.sorted.bam -b -o sample.sorted.unique.bam
        samtools index sample.sorted.unique.bam
step 3: [Predict fragment size]
        python2.7 macs2 predictd -i sample.sorted.unique.bam -f BAM -g genomeSize --outdir sample/
step 4: [Call peak]
        python2.7 macs2 callpeak -t sample.sorted.unique.bam -c control.sorted.unique.bam --extsize fragmentSize --nomodel -f BAM -g genomeSize -n sample -p 1e-5 -B --outdir sample_pe-5/
step 5: [LiftOver]
        liftOver sample.narrowPeak refToGoat.chain sample.goatRef.narrowPeak sample.unmap -minMatch=0.2

Pipeline: ATAC-seq data processing

step 1: [Mapping]
        bowtie2-build ref.fa ref.bowtie2.index
        bowtie2 -x ref.bowtie2.index -X 2000 -1 sample_1.fq.gz -2 sample_2.fq.gz -S sample.sam
        samtools view -bS sample.sam > sample.bam
step 2: [Sort, filter, and remove duplicates and mitochondria from bam]
        samtools sort -@ 8 -o sample.sorted.bam -O bam sample.bam
        samtools view -q 20 sample.sorted.bam -b -o sample.sorted.unique.bam
        java -Xmx50g -Djava.io.tmpdir=tmp -jar picard.jar MarkDuplicates INPUT=sample.sorted.unique.bam OUTPUT=sample.sorted.unique.dedup.bam METRICS_FILE=sample_dedup REMOVE_DUPLICATES=true CREATE_INDEX=true ASSUME_SORTED=true VALIDATION_STRINGENCY=LENIENT MAX_FILE_HANDLES=2000
        samtools view -h sample.sorted.unique.dedup.bam |awk '{if($3 != "Mt_name"){print $0}}' |samtools view -Sb - > sample.sorted.unique.dedup.noMT.bam
step 3: [Call peak]
        python2.7 macs2 callpeak -t sample.sorted.unique.dedup.noMT.bam --nomodel --shift -37 --extsize 73 -f BAM -g genomeSize -n sample -q 0.05 --outdir sample_q0.05/ --broad

Manual

1. Phylogenetic tree

We used whole genome alignments of 51 Ruminantia and one outgroup (killer whale) genomes to construct a high resolution phylogenetic tree. Our results support a sister-group relationship for Antilocapridae and Giraffidae, as well as for Moschidae and Bovidae. Gene flows were observed between Giraffidae and Bovidae, between Giraffidae and Cervidae, and between Giraffidae and Moschidae. Using fossil calibrations, we estimated the emergence of crown Ruminantia at 39.1 - 32.3 Mya (late Oligocene), and the origin of the Pecora at 24.5 - 23.2 Mya (Neocene).

Fig 1. Phylogeny of ruminants.

A. Phylogenetic tree resulted from whole genome alignments of 51 ruminant species using 13 fossil calibrations.

B. Prevalent discordance among single gene trees was observed across different families of ruminant. Left: the 6,406 gene trees were visualized by Densitree; Right: 6,406 simulated gene trees under incomplete lineage sorting (ILS) model, showing indeed extensive ILS exists during the evolution of ruminant families.

2. Families

User can get biological knowledge of six ruminant families from the RGD, including distribution, morphology, and niche.

3. Basic search

3.1 Gene Basic Information and GO/KO Annotation

Users can search with a gene symbol and select one of the above species to get results from three aspects:

(1) "Gene basic information"
    containing Chrom, Start, End, Strand, CDS Length, Exon Number, Transcript ID, Gene Biotype and Transcript
    Product. When selected species is goat, user can click button of "Gbrowse" to jumping the UCSC genome
    browser. Additionally, NCBI Gene link is also provided to help get more information such as
    gene sequence, mRNA sequence and protein sequence.
(2) "Gene Ontology"
    This aspects provides relevant GO ID and GO terms which can be linked to AmiGO 2 database.
    The Gene Ontology is only provided for cattle.
(3) "KEGG pathway"
    Gene pathways linked to KEGG database are also provided for goat, sheep, cattle, and yak. Users can get KO, pathway, motif, protein sequence, nucleic
    acids sequence and etc through this KEGG link.

3.2 Gene Expression Heatmap

Users can input a gene symbol in the search box and choose one species getting a heatmap with different tissues.

(1) Goat [89 samples]:
    Ten tissues at different stages are displayed, especially for rumen at 10 time-points from birth to six month.
(2) Sheep [832 samples]:
    Eighty tissues at different stages covering reproductive, endocrine, immune, integumentary, muscular, embryo, 
    nervous, cardiovascular, respiratory, renal and digestive system, as well as keratin and fat tissues are shown 
    in one heatmap, so user can easily find which tissue has high abundance.
(3) Cattle [163 samples]:
    Forty-one tissues at different stages containing digestive, nervous, immune, reproductive, respiratory, cardiovascular, 
    renal, and muscular system were offered, which may be convenient to study production traits.   
(4) Yak [10 tissues]:
    We offer expression data of 10 tissues (heart, lung, kidney, spleen, brain, overy, muscle, stomach, liver, and ear) for yak, 
    which may be convenient to study adaptation traits.
(4) Roe deer [20 samples]:
    Twenty tissues were collected from a seven-month-old male roe deer. Tissue samples include heart, liver, spleen, 
    lung, kidney, muscle, adipose, cerebrum, testis, esophagus, rumen, reticulum, omasum, abomasum, cecum, ear 
    bone, frontal skin, antler cartilage, antler pedicle and antler velvet. 
(5) Sika deer [20 samples]:
    Twenty tissues were collected from a one-year-old male sika deer. Tissue samples include heart, liver, spleen, 
    lung, kidney, muscle, adipose, cerebrum, cerebellum, testis, esophagus, rumen, reticulum, abomasum, cecum,
    ear bone, frontal skin, antler cartilage, antler pedicle and antler velvet.

4. Tools

4.1 Local UCSC Genome Browser

Users can type a genomic region, or a gene symbol, or a transcript name in the search box and click the "go" button to visualize ruminant data in the UCSC Genome Browser (Gbrowse) . Since comparative genomics data is displayed using goat as reference genome and goat gene annotations are not well-developed, we added four tracks of cross-species gene annotation with "Sheep NCBI Genes", "Cattle NCBI Genes", "Cattle Ensembl Genes", and "hg38 UCSC Genes" (Fig2, Section 1). The Gbrowse also displays an "expreBar" bar graph (Fig2, Section 2) for goat RNA-seq data. A box plot showing the range of expression levels across samples is displayed on the details page when users click the "expreBar" track. Regulatory regions predicted by utilizing 32 ruminant epigenomic and 726 human epigenomic datasets were also displayed (Fig2, Section 3). Moreover, five scales measurements of evolutionary conservation (Fig2, Section 4) and 67-way multiz align (Fig2, Section 5) were also visualized in Gbrowse. Users can find the differences among species by zooming in the alignments to base level. And protein sequences of 67-way alignments can be translated into amino acids. The comparative genomics data can be used as guides to help locate protein-coding genes, identify functional elements, reveal similarities and differences among species, and trace the course of evolution, which will provide more useful information than a single individual genome. Users can click "PDF/PS" item under the "View" menu of navigation bar to generate a high quality image in PostScript or PDF formats.

Figure 2. Main Genome Browser display page on the goat ARS1 assembly, showing gene annotations, expreBar, regulatory regions, eight conservation score tracks, and 67-way multiple alignment track.

4.2 Local UCSC Table Browser

UCSC Table Browser is a powerful tool for retrieving raw data and performing intersections and unions between data in different tracks. For the basic data queries, users can select clade, genome, assembly, group, track, table, regions of interest, output format and output file name to get query results in a tab-delimited text format or compressed format. While for the advanced queries, users can filter and refine queries, intersect query results from different tables and configure the resulting output. The UCSC Table Browser can retrieve and download all data from tracks of goat genome for other analysis.

4.3 Blat

webBlat is a web-based version of BLAT developed by Jim Kent. User can type a DNA or mRNA sequence against the goat ARS1 assembly to return a list of links to all genome positions that share 95% or greater identity with the input sequence. Then the alignment regions can be displayed in the genome browser when users click the "Gbrowse" link.

4.4 Blast

NCBI's wwwBLAST is also available in our database as an online tool, which can be found under the "Tools" menu. User can enter query sequences of DNA against goat genome database getting a group of high-scoring pairwise alignments.