NGS-based-Whole-Genome-Resequencing-Analysis|2Samples
Addendum to the press conference — subject: humpback whale Hope/Timmy
This report analyses the raw DNA data on the humpback whale Hope/Timmy — provided by Dr. med. vet. Kirsten Tönnies together with the initiative #TeamHope.
Who was involved, and where to find them:
- #TeamHopeFacebook TeamHope_Offiziell · TikTok @teamhope_offiziell · Instagram @teamhope_offiziell
- Dr. med. vet. Kirsten TönniesFacebook Tierarztpraxis Dr. Toennies · TikTok @drtoennies · Instagram @drtoennies
- Sergio BambarenFacebook Sergio Bambaren · TikTok @sergiobambaren · Instagram @sergio_bambaren
Data. Transparency. For the truth.
Source: addendum to the press conference, issued by #TeamHope.
Abstract
Summary
2 samples were resequenced. The detailed analysis contents are as follows. Sequencing data evaluation: statistics of sequencing data amount, quality and GC content. Alignment with reference genome: statistics of alignment rate, genome coverage and the depth of genome coverage. Variation identification and annotation: detection and annotation of SNP, InDel, SV and CNV. Variation identification and annotation in coding regions: detection of genes with SNP non-synonymous mutation, InDel mutation, SV mutation and CNV mutation. Gene annotation: annotations on mutated genes by KEGG, GO, COG, NR, SwissProt database.
Overview of results:
132.42Gbp Clean Data was analyzed, and the ratio of Q30 data is 99.14%. The average alignment rate between individual samples and reference genome is 82.66%. The mean coverage depth is 16X. The genome coverage is 93.63% (at least covered by one nucleotide record). The detailed results of identified variations (SNP, Indel) are shown in the report below. Non-synonymous SNPs, InDels, SV and CNV genes were detected between the samples and the reference genome. Gene mutations at DNA level were annotated refer to databases such as KEGG, GO, COG, NR, SwissProt, etc. (for details, the contract shall prevail).
1 General information
Information of reference genome:
Table 1 Sample Information
| Sample ID | fq1 | fq2 |
|---|---|---|
| D1 | D1_1.fq.gz | D1_2.fq.gz |
| T1 | T1_1.fq.gz | T1_2.fq.gz |
Note: Sample ID: Sample ID, which were applied in all steps including library preparation, sequencing and bioinformatics.
Reference Genome:
Table 2 Genome information
| Species | Genome(Mb) | ChrNum | GC(%) | N50(Mb) | Max_Len(Mb) |
|---|---|---|---|---|---|
| Megaptera_novaeangliae | 3146.60 | 1355 | 42.05 | 107.42 | 191.56 |
2 Project process
2.1 Library construction workflow
Library construction workflow follows the standard protocol provided by Illumina Corporation, including sample quality control, library construction, library quality control and library sequencing, etc. The specific flow chart is as follows.
Figure 1 Library construction workflow
After quality control on DNA samples, the libraries were prepared following DNA fragmentation (by sonication), End-repair, adapter ligation on both ends, fragment size selection, PCR enrichment, size selection on PCR products and library quality control. The qualified libraries were sequenced by using Illumina.
2.2 Bioinformatic workflow
The quality of paired-end raw reads were evaluated and the Clean Reads were obtained after filtering and removing low quality reads. By mapping the Clean Reads to the reference genome, genetic variations (like SNPs, InDels .etc) were identified and used to perform DNA-level differential genes mining and differential genes function annotation.
Bioinformatic workflow chart for resequencing is shown as below:
Figure 2 Bioinformatic workflow chart after resequencing
3 Bioinformatic analysis method and result
3.1 Sequencing data quality control
3.1.1 Sequencing data description
Raw image data files obtained from high-throughput sequencing (Illumina sequencing platform)were converted into raw sequenced reads through base calling, which was called as Raw Data or Raw Reads and stored as FASTQ(fq)files, including sequence information and the corresponding sequence quality information. An example of randomly selected real sequencing read is as follows:
Figure 3 Fastq format
Note:Each Read in the FASTQ format file is described by four lines, the first of which begins with “@” followed by Illumina Sequence Identifiers and description (optional part); The second line is the sequence of bases; The third line begins with “+”, followed by the Illumina sequencing identifier (optional part); The fourth line shows the sequencing quality of the corresponding sequence.
Details of the Illumina Sequence Identifiers are as follows:
Table 3 Detailed information table of Illumina Sequence Identifiers
| Identifiers | Detailed Information |
|---|---|
| HWI-7001455 | Unique instrument name |
| 224 | Run ID |
| HA58FADXX | Flowcell ID |
| 2 | Flowcell lane |
| 2216 | Tile number within the flowcell lane |
| 5226 | 'x'-coordinate of the cluster within the tile |
| 55443 | 'y'-coordinate of the cluster within the tile |
| 1 | Member of a pair, 1 or 2 (paired-end or mate-pair reads only) |
| N | Y if the read fails filter (read is bad), N otherwise |
| 0 | 0 when none of the control bits are on, otherwise it is an even number |
| GAAACC | Index sequence |
By using the ASCII value corresponding to each character in the fourth line, the sequencing quality value corresponding to each base in the second line is obtained. If the sequencing error rate is represented by e and the base mass quality value of Illumina is represented by Qphred, then the following relationship exists:
Figure 4 Formula 1
3.1.2 Data quality statistics
Raw Sequenced Reads or Raw Reads were obtained from sequencing, which also have low quality reads and adaptors. To ensure the quality of the data, Raw Reads were filtered to get Clean Reads for subsequent bioinformatic analysis. The main steps of data filtering are as follows : (1) Remove adapter sequences from reads. (2) Remove reads where the number of bases with a quality score Q<15 accounts for more than 40% of the entire read; remove reads with polyG tails; remove reads with sequence complexity < 10%. (3) Remove reads with a length of less than 100 bp. (4) Remove reads with more than 5 N bases in one end
The statistical results of output sequencing data by each sample are shown in the following table:
Table 4 Statistical results of sequencing data by each sample
| _ID | Clean_Reads | Clean_Base | Q20(%) | Q30(%) | GC(%) |
|---|---|---|---|---|---|
| D1 | 355978528 | 53111964460 | 99.82 | 98.98 | 44.35 |
| T1 | 531237660 | 79311686110 | 99.87 | 99.29 | 37.60 |
Note:_ID:The uniform sample ID;Clean_Reads:Number of filtered reads; Clean_Base:number of filtered bases, which equals to the number of Clean_Reads times the sequence length;Q20(%):The percentage of bases with a mass value greater than or equal to 20 in the total number of bases;Q30(%):The percentage of bases with a mass value greater than or equal to 30 in the total number of bases;GC(%):Sample GC content, the percentage of G and C type bases in the total bases.
- 1. PE150 refers to paired-end 150 bp. Why is Clean Bases not a multiple of 150?A: Data quality control procedure to generate clean data is summarized below: 1) Adapter sequence removal: The reads shorter than 100 bp after adapter sequence removal will be removed. 2) N-rich reads removal: N-rich reads indicate those with either one read of the PE reads containing with over 10% of unknown bases (Ns). The paired reads will be removed. 3) Low-quality reads removal: Low-quality reads indicate either one read of the PE reads containing over 50% of bases with quality score lower than Q10. The paired reads will be removed.
3.1.3 Base sequencing quality distribution
The sequencing error rate of each base was obtained by conversion of sequencing Phred score (Qphred) by formula 1, while the Phred value was calculated by a probabilistic model for predicting the occurrence of error in base discrimination during Base Calling. The correlation is shown in the following table:
Table 5 Prediction of error probability for base discrimination
| Phred score | Incorrect base discrimination | Correct base discrimination rate | Q-sorce |
|---|---|---|---|
| 10 | 1/10 | 90% | Q10 |
| 20 | 1/100 | 99% | Q20 |
| 30 | 1/1000 | 99.9% | Q30 |
| 40 | 1/10000 | 99.99% | Q40 |
When applying Illumina sequencing system, the library DNA templates will be firstly fixed on chip. During fixation of DNA template, each DNA molecule will form a cluster and each cluster will become a sequencing site. Very few clusters will have physical overlap. During sequencing process, software will use the first four bases to analyse and distinguish those overlapping sites. Then those sites will be separated to ensure that each sequencing site only have one DNA molecule. The first four cycles are used for cluster calling and establishing more suitable parameter setting based on current flowcell. Thus, the error rates of first few bases at 5' end of a sequence are relatively higher. In addition, the sequencing error rate will raise as the length of sequenced read increases, which is caused by the consumption of chemical reagents in the sequencing process. Therefore, the quality values of the first four and the last a dozen or so bases are lower than those of intermediate sequencing bases. Still, the quality values would be higher than Q30%. According to the relationship between the quality value and the error rate, we converted the quality value into the error rate and plotted the error rate distribution diagram as follows:
Figure 5 Sample base error rate distribution
Note:x-axis represents base position of reads, y-axis represents the error rate of single base. The first 150bp shows the error rate distribution of the first end of paired-end sequences, the second 150bp shows the error rate distribution of another end.)
3.1.4 Base composition distribution
Base composition distribution check is used to detect whether there is AT or GC separation caused by sequencing or library construction, which may affect subsequent analysis. Given the randomness of the sequence and base complementarity, distribution of base content should be approximately even and the contents of G/C and A/T are theoretically equal on each sequencing circle, and also equal to the GC and AT contents of the genome respectively. Similarly, due to cluster calling, the AT and GC content on the first few bases are highly fluctuant. While in other regions, the AT and GC content are evenly ditributed without seperation, as shown in the figure below:
Figure 6 Proportion distribution of each base of the sample
Note:x-axis represents base position of reads, y-axis represents the proportion of bases; Different colors represent different base types, green for base G, blue for base C, red for base A, purple for base T, and gray for base N, which is not recognized by sequencing. The first 150bp shows the base distribution of the first end of paired-end sequences. The second 150bp shows the base distribution of another end. Each cycle represents each base sequenced. For example, the first cycle represents the distribution of A, T, G, C and N of all sequencing reads at the first base. This figure shows that the AT and CG bases are basically not separated and the curve is relatively flat: So the sequencing result is normal.
3.2 Statistics of alignment to reference genome
Sequencing reads obtained by resequencing need to be mapped against the reference genome for subsequent variation analysis. The bwa[1] software is mainly used to align the short sequences obtained by second-generation high-throughput sequencing (Sequencing platforms such as Illumina) to reference genome. By aligning Clean Reads to the reference genome, the sequencing depth, genome coverage and other information of each sample were collected, and the variation was detected.
3.2.1 Statistics of alignment results
The statistics of reference genome alignment including Mapping ratio、Sequencing depth and coverage were performed by samtools flagstat/depth.
Reads were aligned to reference genome to obtain information. Alignment rate: the proportion of Clean Reads aligned to the reference genome to the total Clean Reads. If the appropriate reference genome is selected and there is no contamination in experimental process, the alignment rate of sequencing reads will be higher than 70%. In addition, the alignment efficiency is related to the genetic distance between sequenced species and the reference genome, the assembly quality of reference genome and the sequencing reads quality. For instance, closer genetic distance, higher integrity of reference genome, and higher quality of the sequencing reads, will help map more reads to the reference genome and improve the alignment rate.
Sample alignment result is shown in the following table:
Table 6 Statistics of alignment results
| ID | Total_reads | Mapped(%) | Properly_mapped(%) | Ave_depth | Cov_ratio_1X(%) | Cov_ratio_5X(%) | Cov_ratio_10X(%) |
|---|---|---|---|---|---|---|---|
| D1 | 355978528 | 79.50 | 73.61 | 12 | 90.77 | 66.48 | 46.17 |
| T1 | 531237660 | 85.81 | 83.67 | 20 | 96.5 | 85.54 | 70.52 |
Note:ID:The uniform sample ID;Total_reads:number of Clean Reads;Mapped(%): Percentage of Clean Reads mapped to the reference genome;Properly_mapped(%):Percentage of paired-end Clean Reads mapped to the reference genome with proper distance refers to sample fragment length in between.
Coverage and sequencing depth on each chromosome were visualized as shown below. In order to guarantee an optimized appearance, a maximum of 20 chromosomes/scaffolds will be displayed in the report by default. An even distribution of sequencing reads on reference genome indicates a good randomness of sequencing data.
Figure 7 Coverage and sequencing depth
Note: X-axis: Position on chromosomes. Y-axis: Log2 of sequencing depth at corresponding position.
The above figures show an overall evenly distributed coverage on the reference genome, which indicated good randomness of sequencing data. The regions with uneven colour may be repetitive sequences, or bias caused by PCR. If there are regions without any reads aligned, it may be because of the gaps in reference genome.
3.2.2 Statistical distribution of insertion size
By detecting the start and end position of the paired-end sequences which are mapped against the reference genome, the actual size of sequenced fragments from sample DNA can be obtained, i.e. the inserted DNA size, which is an important parameter in bioinformatics. The distribution of insertion size generally conforms to the normal distribution, it has only one single peak, the insertion size distribution map shows the length distribution of insertions of each sample. The analysis of insertion size of sequencing data of each sample is realized by CollectInsertSizeMetric.jar from picard software tool pack.
Figure 8 Insertion distribution map
Note:x-axis represents the length of insertions, y-axis represents number of its corresponding reads.
As the above figure shows, the insertion length distribution fits the normal distribution, indicating that there is no anomaly in the library construction of sequencing data.
3.2.3 Statistics of depth distribution
After mapping Reads to the reference genome, the base coverage on the reference genome can be calculated. The percentage of bases on the reference genome covered by reads is called genome coverage; The number of reads covering the same bases is called the coverage depth. Genome coverage can reflect the integrity of Variation identification. The more regions covered, the more mutated sites that can be detected. The coverage was mainly affected by the sequencing depth and the genetic distance between the sample's genome and the reference genome. The coverage depth of the genome will affect the accuracy of Variation identification. Higher coverage depth guarantees higher accuracy of Variation identification (non-repetitive sequence region). In addition, if the distribution of coverage depth per base on the genome is relatively uniform, it also indicates the randomness of sequencing. The distribution curves of coverage depth per base and coverage degree are shown in the following figure:
Figure 9 Sample depth distribution
Note:The above figure reflects sequencing depth distribution, the x-axis represents sequencing depth, the left y-axis represents the percentage of bases corresponding to this depth, which corresponds to the red curve. The right y-axis represents the percentage of cumulative bases, which corresponds to the blue curve.
- 1. What is the normal range of mapping ratio [Mapped(%)]?A.With a good quality of data and a proper reference genome, the mapping ratio should be no smaller than 70%.
- 2. What are the potential reasons for a low mapping ratio?A: 1)The quality of reference genome is poor. 2) The species of reference genome is not close to that of the materials. 3) Material has be processed for experiments which caused changes in genome DNA. 4) The samples are contaminated by other species.
- 3. How to deal with data with low mapping ratio?A: 1) Re-check the quality of reference genome. If the assembly used in currently analysis is of low-quality, we recommend to search for another assembly version. 2) Clarify the genetic relationship between sequencing material and reference genome. 3) Check if there are contamination from other species using NT database. The contaminated reads can be removed accordingly. If valid sequencing depth is not sufficient after data filtration, we recommend to re-sequence the same library to add on more data.
- 4. What are the potential reasons for a low coverage [Cov_ratio_1X(%)]?A: Normally, a dataset with an average of 10X sequencing depth can achieve 1X coverage of over 90%. Potential reasons of coverage lower than that can be: 1) Low mapping ratio, i.e. the sequencing depth of mapped reads is lower than expectation. 2) Large difference in reference genome and sequencing material.
3.3 Variation identification and annotation
3.3.1 Tools and methods for Variation identification
Detecting SNP(Single Nucleotide Polymorphism)and small InDel(small Insertion and Deletion)were mainly conducted by GATK[2] software package. According to the mapping result of Clean Reads to the reference genome, samtools(v1.9)[3] was used to filter abundant reads(MarkDuplicates)and ensure the accuracy of testing result. Then the HaplotypeCaller(Local haplotype assembly)algorithm from GATK was applied to detect SNP and InDel variations.It firstly generated gVCF files from each sample, then performed population joint-genotype . Finally, mutated sites set was collected after filtering.
Please find the website for the detailed description of vcf file:https://github.com/samtools/hts-specs
3.3.2 Quality control of Variation identification result
Variants were filtered strictly to ensure its credibility, the main filtering parameters are as follows:
(1)Use subprogram vcfutils.pl(varFilter -w 5 -W 10)from bcftools to filter out SNPs in the vicinity of INDEL within 5bp and SNPs in the adjacent INDEL within 10bp;
(2)clusterSize 2 clusterWindowSize 5. This parameter indicates that the number of variants in the 5bp window should not exceed 2;
(3)QUAL < 30. This parameter indicates that variants whose quality score in Phred format lower than 30 will be filtered out;
(4)QD < 2.0 The parameter QD represents the ratio of the variant's quality score divided by its coverage depth. The variant's coverage depth is the sum of the samples at this site containing variant bases. Variants with QD value less than 2.0 are filtered out;
(5)MQ < 40. MQ represents the Root Mean Square of alignment quality value of the reads aligned to this site. Values with MQ less than 40 are filtered out;
(6)FS > 60.0. FS value is converted from the p-value of Fisher's test. It describes whether there is an significant positive or negative chain specificity to the reads that only contains reference genotype or mutated genotype. In other words, FS should be close to zero if there is no chain-specific alignment. Thus, values with FS above 60 are filtered out;
(7)Other variation filtering parameters are treated by following the default settings of GATK.
- 1. How can I validate the SNPs identified in the data?A: The most commonly used method for SNP validation is PCR and KASP.
3.3.3 High quality variation demonstration
SNP mutations can be divided into transitions and transversions. Mutations between the same type of bases are named transitions, such as those between purines and purines and between pyrimidines and pyrimidines. Mutations between different types of bases are named transversions, i.e. those between purines and pyrimidines. Generally transition is more likely to occur comparing to transversion. Therefore, the ratio of transition/transversion (Ti/Tv) is normally greater than 1. This value is dependent on specific species. For diploid or polyploid species, if a certain SNP site on homologous chromosomes is the same base, this SNP site is called homozygous SNP site; If the SNP site on homologous chromosomes contains different types of bases, the SNP site is called heterozygous SNP site. The higher the number of homozygous SNP sites, the greater the difference between the sample and the reference genome. The higher the number of heterozygous SNP sites, the higher the heterozygosity of the sample. The SNP detection results between the sample and the reference genome were shown below.
Table 7 Statistics on SNP identification
| ID | SNPnumber | Transition | Transversion | Ti/Tv | Heterozygosity | Homozygosity | Het-ratio |
|---|---|---|---|---|---|---|---|
| D1 | 3707264 | 2535550 | 1171714 | 2.16 | 2176650 | 1530614 | 58.71% |
| T1 | 4309415 | 2921983 | 1387432 | 2.1 | 2623529 | 1685886 | 60.87% |
| Total | 4895867 | 3300981 | 1594886 | 2.06 | - | - | - |
Note: ID: Sample ID. SNPnumber: Counts of identified SNPs. Transition: Counts of transition SNPs. Transversion: Counts of transversion SNPs. Ti/Tv: The ratio of transition SNPs to transversion SNPs. Heterozygosity: Counts of heterozygous SNPs. Homozygosity: Counts of homozygous SNPs. Het-ratio: The percentage of heterozygous SNPs.
To ensure the credibility of SNP, collect the cumulative distribution of number of reads supporting the detected SNP and the distance between adjacent SNPs.
Figure 10 SNP quality distribution map
Note:The left figure is the cumulative distribution map of reads supporting the SNP, The right figure is the cumulative map of the distance between adjacent SNPs.
3.3.4 SNP detection between samples (multiple samples)
According to the alignment result between samples and the reference genome, the variation sites of all samples are summarized. The format of the SNP list file of all samples is as follows:
Table 8 SNP list of samples
| #Chr | Pos | Ref | D1 | T1 |
|---|---|---|---|---|
| 1 | 966 | T | T | Y |
| 1 | 1247 | C | C | S |
| 1 | 1749 | G | G | R |
| 1 | 2589 | G | R | G |
| 1 | 2601 | G | K | G |
Please find the detailed data here:Up_load/03.SNP/04.DEG/Megaptera_novaeangliae.raw.filter.snp.DEG.snp
SNP genotypes are encoded using standard nucleotide symbols, as shown in the symbol table below:
Table 9 SNP genotypes code
| Nucleotide Symbols | Meaning | Nucleotide Symbols | Meaning |
|---|---|---|---|
| A | Adenosine | M | A C (aMino group) |
| C | Cytosine | S | G C (Strong interaction) |
| G | Guanine | W | A T (Weak interaction) |
| T | Thymidine | B | G T C (not A) (B comes after A) |
| U | Uracil | D | G A T (not C) (D comes after C) |
| R | G A (puRine) | H | A C T (not G) (H comes after G) |
| Y | T C (pYrimidine) | V | G C A (not T, not U) (V comes after U) |
| K | G T (Ketone) | N | A G C T (aNy) |
Genome-wide SNP mutations can be divided into six categories. Using T:A>C:G as an example, this type of SNP mutation has T>C and A>G. Since the sequencing data can be aligned to the positive chain or the negative chain of the reference genome, when T>C type mutation occurs at the positive chain of the reference genome, A>G type mutation will occur at the same position of the genome negative chain. So T>C and A>G are classified as one class.
Figure 11 Distribution of SNP types
Note:y-axis represents number of SNPs, x-axis represents SNP mutation type.
The statistical results of SNPs among samples are shown in the following figure:
Figure 12 Statistical Venn diagam of SNP among samples
- 1. After verifying the variant sites in the DEG results, it was found that they were not unique to the sample.A: Whether a variant site is unique or not is determined based on complete genotype matching. For example, if sample A has a genotype of AA and sample B has AT, this is considered a differential site. It is recommended to verify homozygous and differential sites, i.e. sites whose genotypes are A, T, C, and G after degeneracy conversion.
3.3.5 SNP annotation
SnpEff[4] is a software used to annotate variations (SNP, Small InDel) and predict the effect of them. According to the position of the mutated locus and information on the reference genome , the region where the mutated locus occurs (intergene region, gene region or CDS region, etc.) and the impact of the mutation (synonymous mutation, non-synonymous mutation, etc.) can be obtained. The software can use VCF format files as input and output. The output result will add the following fields to the INFO column of the VCF file:EFF= Effect ( Effect_Impact | Functional_Class | Codon_Change | Amino_Acid_Change| Amino_Acid_Length | Gene_Name | Transcript_BioType | Gene_Coding | Transcript_ID | Exon_Rank | Genotype_Number [ | ERRORS | WARNINGS ] )
SNP classification of each sample were shown in the following table.
Table 10 Statistics on SNP annotation
| #Region | Type | All | D1 | T1 |
|---|---|---|---|---|
| -- | INTERGENIC | 2551867 | 1870731 | 2252936 |
| -- | INTRAGENIC | 14849 | 11765 | 13227 |
| -- | INTRON | 2007478 | 1562142 | 1772391 |
| -- | UPSTREAM | 119359 | 101508 | 100672 |
| -- | DOWNSTREAM | 114138 | 95201 | 98863 |
| -- | UTR_5_PRIME | 6036 | 5521 | 3355 |
| -- | UTR_3_PRIME | 28230 | 22694 | 24435 |
| -- | SPLICE_SITE_ACCEPTOR | 89 | 66 | 73 |
| -- | SPLICE_SITE_DONOR | 141 | 120 | 111 |
| -- | SPLICE_SITE_REGION | 3028 | 2562 | 2395 |
| -- | START_GAINED | 396 | 343 | 320 |
| CDS | START_LOST | 41 | 36 | 30 |
| CDS | SYNONYMOUS_START | 0 | 0 | 0 |
| CDS | NON_SYNONYMOUS_START | 0 | 0 | 0 |
| CDS | SYNONYMOUS_CODING | 30183 | 18239 | 25068 |
| CDS | NON_SYNONYMOUS_CODING | 19725 | 16128 | 15293 |
| CDS | SYNONYMOUS_STOP | 22 | 18 | 14 |
| CDS | STOP_GAINED | 255 | 163 | 207 |
| CDS | STOP_LOST | 30 | 27 | 25 |
| -- | Other | 0 | 0 | 0 |
Distribution of SNPs is visualized in the pie charts below.
Figure 13 Pie chart of SNP annotation
Terminologies in SNP annotation were summarized in the following table.
Table 11 Functional regions
| Function or Type | Meaning |
|---|---|
| INTERGENIC | Intergenic region |
| INTRAGENIC | Intragenic(no transcript information) |
| INTRON | Intron |
| UPSTREAM | Upstream region of gene (within 5K) |
| DOWNSTREAM | Downstream region of gene (within 5K) |
| UTR_5_PRIME | 5’UTR region |
| UTR_3_PRIME | 3’UTR region |
| SPLICE_SITE_ACCEPTOR | Splicing acceptor variation (within the first 2bp of exon) |
| SPLICE_SITE_DONOR | Splicing donor variation (within the last 2bp of exon) |
| SPLICE_SITE_REGION | Splicing region variation (1-3 bases variation in exon or 3-8 bases variation in intron) |
| START_GAINED | Start codon gain (non-coding regions) |
| START_LOST | Start codon lost |
| SYNONYMOUS_START | Synonymous start codon variation |
| NON_SYNONYMOUS_START | Nonsynonymous start codon variation |
| SYNONYMOUS_CODING | Synonymous coding variation |
| NON_SYNONYMOUS_CODING | Nonsynonymous coding variation |
| SYNONYMOUS_STOP | Synonymous stop codon variation |
| STOP_GAINED | Stop codon gain |
| STOP_LOST | Stop codon lost |
| Other | Due to the incomplete genetic information in the GFF file, accurate information cannot be obtained |
3.3.6 Small InDel identification between sample and the reference genome
According to the mapping results of sample Clean Reads on the reference genome, the insertion and deletion of small fragments (Small InDel) were detected. GATK was used to detect InDel of the sample. The number of Small InDel variations is generally less than those of SNP variations, which also reflects the difference between the sample and the reference genome. InDels in the coding region will cause code shift mutation, leading to changes in gene function. Please find the table below
Table 12 Statistics of InDels in whole genome and the coding region
| #Sample | CDS-Insertion | CDS-Deletion | CDS-Homo | CDS-Het | CDS-Total | Genome-Insertion | Genome-Deletion | Genome-Homo | Genome-Het | Genome-Total |
|---|---|---|---|---|---|---|---|---|---|---|
| D1 | 763 | 1003 | 370 | 1396 | 1766 | 233920 | 293371 | 198505 | 328786 | 527291 |
| T1 | 1065 | 3709 | 485 | 4289 | 4774 | 340075 | 537938 | 236386 | 641627 | 878013 |
| #Total | 1623 | 4437 | -- | -- | 6060 | 407994 | 635744 | -- | -- | 1043738 |
Note:CDS:coding region InDel statistics;Genome:Genome-wide InDel statistics;Insertion:Number of detected insertions;Deletion:Number of detected deletions;Het:Number of heterozygous InDel;Homo:Number of homozygous InDel;Total:Total number of detected InDel (Duplicates removed).
Statistics of the InDel length in the CDS region and the genome-wide region, the length distribution is shown in the figure below(in order to guarantee an optimized appearance,only top 50 samples will be displayed in the figure by default):
Figure 14 Distribution map of InDel of whole genome and coding region
Note:Y-axis represents the length of InDel (within 10bp). Positive value means Insertion and negative value means Deletion. X-axis represents the corresponding number.
3.3.7 Small InDel identification between samples
Part of the Small InDel identification results between samples and reference genome are shown in the following table.
Table 13 Statistics of Small InDel sequencing data of samples
| #Chr | Pos | Ref | Alt | D1 | Depth | AlleDp | T1 | Depth | AlleDp | Effect | Codon_change | Gene_id | COG_class | COG_class_annotation | GO_annotation | KEGG_annotation | KOG_class | KOG_class_annotation | Pfam_annotation | Swissprot_annotation | TrEMBL_annotation | eggNOG_class | eggNOG_class_annotation | nr_annotation |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 622 | GA | G | GA,G | 15 | 11,4 | GA,G | 18 | 7,11 | INTERGENIC | - | |||||||||||||
| 1 | 942 | C | CAA | C,C | 9 | 9,0 | C,CAA | 13 | 10,3 | INTERGENIC | - | |||||||||||||
| 1 | 1633 | CTG | C | CTG,C | 18 | 15,3 | CTG,C | 33 | 23,10 | INTERGENIC | - | |||||||||||||
| 1 | 2534 | GA | G | GA,G | 31 | 26,5 | GA,G | 23 | 21,2 | INTERGENIC | - | |||||||||||||
| 1 | 3826 | TC | T | TC,TC | 19 | 19,0 | TC,T | 7 | 5,2 | INTERGENIC | - |
Note: #Chr: Chromosome ID; Pos: Position in reference genome; Ref: Reference base(s); Alt: Alternate base(s) of InDel; The following three columns are genotype, sequencing depth of each genotype and sequencing depth of the site. Effect: InDel type provided by SnpEff. Codon_change: Changed codon due to frameshift; Gene_id: Gene affected by the InDel. The following columns are annotation of this gene in each database.
The statistical result of INDEL between samples is shown in the following figure:
Figure 15 Statistical Venn diagam of INDELs among samples
3.3.8 Small InDel annotation
According to the location information of detected Small InDel and information such as reference gene and CDS location, etc. (generally can be found in gff files), one can annotate whether the InDel occurs in the intergene region, gene region or CDS region, or whether it is a code shift mutation. The annotation of Small InDels is conducted by SnpEff software. InDels with code shift mutation may result in changes of gene function. The detailed annotation result is shown in the figure below:
Table 14 Statistics on InDel annotation
| #Region | Type | All | D1 | T1 |
|---|---|---|---|---|
| -- | INTERGENIC | 507976 | 245540 | 431144 |
| -- | INTRAGENIC | 3369 | 1738 | 2815 |
| -- | INTRON | 455836 | 238187 | 382409 |
| -- | UPSTREAM | 30489 | 18014 | 24000 |
| -- | DOWNSTREAM | 29099 | 16334 | 24206 |
| -- | UTR_5_PRIME | 1404 | 969 | 728 |
| -- | UTR_3_PRIME | 8354 | 4243 | 7017 |
| -- | SPLICE_SITE_ACCEPTOR | 74 | 26 | 62 |
| -- | SPLICE_SITE_DONOR | 114 | 47 | 86 |
| -- | SPLICE_SITE_REGION | 963 | 427 | 772 |
| CDS | START_LOST | 2 | 2 | 1 |
| CDS | FRAME_SHIFT | 4786 | 1024 | 4032 |
| CDS | CODON_DELETION | 328 | 223 | 172 |
| CDS | EXON_DELETED | 1 | 1 | 1 |
| CDS | CODON_INSERTION | 340 | 221 | 191 |
| CDS | CODON_CHANGE_PLUS_CODON_DELETION | 249 | 147 | 144 |
| CDS | CODON_CHANGE_PLUS_CODON_INSERTION | 162 | 87 | 96 |
| CDS | STOP_GAINED | 184 | 55 | 131 |
| CDS | STOP_LOST | 8 | 6 | 6 |
| -- | Other | 0 | 0 | 0 |
Distribution of InDels is visualized in the pie charts below.
Figure 16 Pie chart of InDel annotation
Terminologies in InDel annotation were summarized in the following table.
Table 15 Terminologies in SNP annotation were summarized in the following table.
| Function or Type | Meaning |
|---|---|
| INTERGENIC | Intergenic region |
| INTRAGENIC | Intragenic(no transcript information) |
| INTRON | Intron |
| UPSTREAM | Upstream region of gene (within 5K) |
| DOWNSTREAM | Downstream region of gene (within 5K) |
| UTR_5_PRIME | 5’UTR region |
| UTR_3_PRIME | 3’UTR region |
| SPLICE_SITE_ACCEPTOR | Splicing acceptor variation (within the first 2bp of exon) |
| SPLICE_SITE_DONOR | Splicing donor variation (within the last 2bp of exon) |
| SPLICE_SITE_REGION | Splicing region variation (1-3 bases variation in exon or 3-8 bases variation in intron) |
| START_LOST | Start codon lost |
| FRAME_SHIFT | Frame shift variation (insert or delete multiple of non-3) |
| CODON_DELETION | Codon deletion (integer multiple of 3) |
| EXON_DELETED | Exon deleted |
| CODON_INSERTION | Codon insertionn (integer multiple of 3) |
| CODON_CHANGE_PLUS_CODON_DELETION | Codon change + codon deletion |
| CODON_CHANGE_PLUS_CODON_INSERTION | Codon change + codon insertion |
| STOP_GAINED | Stop codon gain |
| STOP_LOST | Stop codon lost |
| Other | Due to the incomplete genetic information in the GFF file, accurate information cannot be obtained |
3.3.9 SV identification
Genomic structural variation (SV) refers to the variation of chromosomal structure of species, such as large fragment insertion, deletion, inversion and translocation. Manta is commonly used to detect structural variations. It firstly obtained the inserted DNA fragment size and variance of the sequencing library from the alignment against the reference genome. Then it identified possible structural variations by looking for abnormal alignment results between the sequence and the reference genome (deviation of the insertion fragment, inconsistent alignment direction, etc.). The SV result file contains the title line and the data line, as shown below:
Manta[5] software detect insertion (INS), Deletion (DEL), Inversion (INV), Chromosomal Translocation (TRA) by checking the alignment between the paired-end reads and the reference genome and the actual inserted DNA size. Statistics of SV numbers and for each type of samples is shown in the following table:
Table 16 Demo result of SV calling
| #Sample | Chr1 | Pos1 | Orientation1 | Chr2 | Pos2 | Orientation2 | Type | Size | Score | SV_loc | Gene |
|---|---|---|---|---|---|---|---|---|---|---|---|
| D1 | 1 | 506871 | - | 1 | 506956 | - | INV | 85 | , | gene | gene-AAG666_011756 |
| D1 | 1 | 654404 | - | 1 | 654488 | - | INV | 84 | , | intergenic | NA |
| D1 | 1 | 662921 | - | 1 | 663035 | - | INV | 114 | , | gene | gene-AAG666_014865 |
| D1 | 1 | 1128196 | - | 1 | 1128254 | - | DEL | 58 | , | gene | gene-AAG666_018440 |
| D1 | 1 | 1164617 | - | 1 | 1164779 | - | DEL | 162 | , | intergenic | NA |
Note: Chr1: Chromosome ID, where the first breakpoint of SV is located. Pos1: Position of the first breakpoint; Orientation1: Number of reads mapped to the first breakpoint (+ * sense reads, - * ani-sense reads); Chr2: Chromosome ID, where the second breakpoint of SV is located. Pos2: Position of the second breakpoint; Orientation2: Number of reads mapped to the second breakpoint (+ * sense reads, - * ani-sense reads); Type: SV type; Size: Length of SV; Score: Confidence value; SV_loc: SV position; Gene: Annotated genes.
Statistics of identified SVs in 2 samples were summarized below.
Table 17 Statistical table SV numbers
| ID | SV_total | INS | DEL | INV | DUP | TRA |
|---|---|---|---|---|---|---|
| D1 | 7252 | 285 | 2203 | 3450 | 310 | 1004 |
| T1 | 7131 | 566 | 4988 | 134 | 241 | 1202 |
Note:ID:The uniform sample ID;SV: Total number of structural variations;INS:number of insertion variations;DEL:number of deletion variations;INV:number of invertion variations;DUP:number of duplication variations;TRA:number of chromosomal translocation variations.
3.4 CNV identification
Use FREEC[6] to detect CNV by the depth distribution of sample sequencing reads on the reference genome. Then draw distribution of Copy Number Gain and Loss on the reference genome. Part of FREEC result file of sample D1:
Table 18 FREEC result file
| chr | start | end | predicted copy number | type of alteration |
|---|---|---|---|---|
| 1 | 62850000 | 62900000 | 5 | gain |
| 1 | 69900000 | 69950000 | 6 | gain |
| 1 | 122850000 | 123900000 | 1 | loss |
| 1 | 131850000 | 131900000 | 5 | gain |
| 1 | 135350000 | 135550000 | 1 | loss |
Note:chr:chromosome;strat:the start position of CNV ;end:the end position of CNV;predicted copy number:predicted copy number;type of alteration:CNV type.
3.5 Distribution of different types of variation on genome
The distribution of different types of variation was shown in the circos plot, which was made by circos software. Website of the software:http://circos.ca/. The distribution of each type of variation on the chromosome of the all samples is shown in the figure below:
Figure 17 The distribution of different types of variation on genome
Note:From the peripheral to the central layers are: chromosome coordinates,Gene density distribution,SNP density distribution, InDel density distribution.
3.6 Gene mutation analysis on DNA level
3.6.1 Gene mutation mining on DNA level
Mutations occurring in the CDS region may cause change of gene function. By targeting non-synonymous mutation SNP, InDel and SV genes occurring in the CDS region, genes with possible functional differences between the samples and the reference genome can be detected. The variation between each sample and the reference genome is shown in the following table:
Table 19 Taxonomic statistics of differential gene mutations
| _ID | Genes with Non-synonymous SNP | Genes with InDel | Genes with SV |
|---|---|---|---|
| D1 | 8731 | 1851 | 2159 |
| T1 | 8665 | 4200 | 1786 |
Note:_ID:The uniform sample ID;Genes with Non-synonymous:The number of genes with non-synonymous mutations. The case that one gene with multiple non-synonymous mutations were not counted repeatedly;Genes with InDel:Number of genes with Small InDel;Genes with SV:Number of genes with SV.
COG database:It is a database for orthologous classification of gene products. Each COG protein is assumed to be from an ancestor protein and is divided into orthologs and paralogs. Orthologs are proteins from different species that have evolved from vertical families and typically retain the same functions as the original proteins. Paralogs are proteins derived from gene replication in a given species that may evolve new functions related to the original.
KEGG database:KEGG integrates the current knowledge of chemical compounds, reactions and molecular networks in biochemistry. It is a database for the systematic analysis of the metabolic pathways and functions of gene products in cells, which is beneficial to integrate gene and expression information as a whole network study.
The annotation list of gene mutation is shown in the table below:
Table 20 Annotation list of mutated genes of sample D1
| #GeneID | COG_class | COG_class_annotation | GO_annotation | KEGG_annotation | KOG_class | KOG_class_annotation | Pfam_annotation | Swissprot_annotation | TrEMBL_annotation | eggNOG_class | eggNOG_class_annotation | nr_annotation |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| gene-AAG666_017024 | -- | -- | Molecular Function: DNA-binding transcription factor activity, RNA polymerase II-specific (GO:000098... | ZBTB40; zinc finger and BTB domain-containing protein 40 | [R] | General function prediction only | Zinc finger, C2H2 type | Telomere zinc finger-associated protein OS=Pongo abelii OX=9601 GN=ZBTB48 PE=2 SV=1 | zinc finger and BTB domain-containing protein 40 OS=Balaenoptera acutorostrata scammoni OX=310752 GN=ZBTB40 PE=4 SV=1 | S | Function unknown | zinc finger and BTB domain-containing protein 40 [Balaenoptera acutorostrata scammoni] |
| gene-AAG666_018102 | -- | -- | Molecular Function: GTPase activator activity (GO:0005096);; Biological Process: signal transduction (GO:0007165); | ARHGAP20; Rho GTPase-activating protein 20 | [T] | Signal transduction mechanisms | RhoGAP domain | -- | rho GTPase-activating protein 20 isoform X1 OS=Balaenoptera acutorostrata scammoni OX=310752 GN=ARHGAP20 PE=4 SV=1 | T | Signal transduction mechanisms | rho GTPase-activating protein 20 isoform X1 [Balaenoptera acutorostrata scammoni] |
| gene-AAG666_002100 | -- | -- | Molecular Function: DNA-binding transcription factor activity, RNA polymerase II-specific (GO:000098... | -- | [R] | General function prediction only | Zinc-finger double domain | -- | zinc finger protein 696 isoform X1 OS=Physeter macrocephalus OX=9755 GN=ZNF696 PE=4 SV=1 | K | Transcription | zinc finger protein 696 [Phocoena sinus] |
| gene-AAG666_025815 | -- | -- | Biological Process: DNA damage checkpoint signaling (GO:0000077);; Biological Process: DNA repair (GO:0006281);; Cellular Component: checkpoint clamp complex (GO:0030896); | RAD9B; cell cycle checkpoint control protein RAD9B | [DL] | -- | Rad9 | Cell cycle checkpoint control protein RAD9B OS=Bos taurus OX=9913 GN=RAD9B PE=2 SV=1 | Cell cycle checkpoint control protein OS=Balaenoptera acutorostrata scammoni OX=310752 GN=RAD9B PE=3 SV=1 | -- | -- | cell cycle checkpoint control protein RAD9B isoform X1 [Balaenoptera acutorostrata scammoni] |
| gene-AAG666_002612 | [R] | General function prediction only | Molecular Function: guanyl-nucleotide exchange factor activity (GO:0005085);; Cellular Component: Go... | ARFGEF, BIG; brefeldin A-inhibited guanine nucleotide-exchange protein | [U] | Intracellular trafficking, secretion, and vesicular transport | Sec7 domain | Brefeldin A-inhibited guanine nucleotide-exchange protein 1 OS=Bos taurus OX=9913 GN=ARFGEF1 PE=1 SV=1 | brefeldin A-inhibited guanine nucleotide-exchange protein 1 isoform X1 OS=Balaenoptera musculus OX=9771 GN=ARFGEF1 PE=4 SV=1 | U | Intracellular trafficking, secretion, and vesicular transport | brefeldin A-inhibited guanine nucleotide-exchange protein 1 isoform X1 [Balaenoptera musculus] |
Statistical result of GO classification on gene mutation is shown in below figure:
Figure 18 GO annotation and clustering of mutated genes of sample D1
Note:x-axis represents each classification of GO. The left y-axis represents the percentage of genes, the right shows number of genes
Statistical result of COG classification of gene mutation is shown in below figure:
Figure 19 COG classification map of gene mutation of sample D1
Note:x-axis represents each classification of COG. Y-axis represents number of genes in different functional classes. The proportion of genes reflects the metabolic or physiological bias in the corresponding period and environment, which can be explained scientifically in combination with the distribution of research objects in each functional class.
Result of metabolic pathways of gene mutation is shown in the figure below:
Figure 20 Metabolic pathways map of gene mutation of sample D1
Note:Numbers in box stand for enzyme ID numbers, suggesting that the corresponding gene is associated with the enzyme. The whole pathway is formed by a variety of different enzymes through complex biochemical reactions. All mutated genes related to this pathway are shown in red boxes. Researchers can focus on the genes of related metabolic pathways according to their research objects, and explain the origin of the corresponding genes through metabolism.
4 Appendix
4.1 Result Directory Readme
Result Directory Readme can be checked in Appendix 1. Result Directory Readme
4.2 Software and Database
Result Directory Readme can be checked in Appendix 2. Parameters of Software and Database
4.3 Material and Method
For the convenience of writing articles, the material and method in Chinese and English involved in the data analysis can be checked in Appendix 3. Material and Method
4.4 Glossary
The technical terms mentioned in the report can be checked in Appendix 4. Glossary
4.5 Acknowledgments and Citations
4.6 Description for viewing result files
1)File decompression: Linux system decompression: zip pack decompression command 'unzip file.zip'; tar pack decompression command 'tar-xvf file.tar'; tar.gz pack decompression command 'tar-xzvf file.tar.gz'; Windows system can use WINRAR or 7-zip for decompression.
2)View text files: Linux or Unix users can use 'more', 'less' and other commands to view; Windows users can use a text editor or Wordpad to open and edit text files, such as gedit or editplus.
3)Image viewing and editing: some image files may be included in the data. Image files generally have suffix such as .png, .pdf, tiff, svg, etc. For image files, Windows users can use the image browser to open them, Linux/Unix users can use display command to open them. If image need to be edited, it is recommended to use the AI editor. Detailed operation: firstly, install Adobe Illustrator (Adobe Illustrator CS6) software on Windows. Double-click .exe file to open it; Secondly, use the shortcut keys 'Ctrl+o' or select 'file-open' to find the PDF file and open it. If it’s single page, the file can be edited directly; if it’s multipage PDF, select the page that need to be edited and then continue to select the left toolbar for editing; Thirdly, save the file. To save the original multipage PDF, simply use the 'save' command, or use shortcut keys 'Ctrl+s'. If you want to save the current editing page as single page PDF, use 'save as' command, or the shortcut keys 'Shift+Ctrl+S'.
4)Open table: the table under Linux is Tab-segmented text, which can be viewed directly with the less command (less-sn *.xls) or opened by excel or openoffice or other office software.
(note: opening a large file can cause Windows crash. It is recommended to use a more powerful computer or a Unix/Linux system that is better suited for handling large amounts of data.)
4.7 FAQ
1. Why some links in the report cannot be opened?
Answer: Links in the report link to all analysis results and original data. After receiving these results, you can open the links to view a analysis summary.
2. Where can I find high-resolution or vector graphics of the images in the report?
Answer: Report images are always in PNG or JPG format with two conversions to control the scale of the report. Corresponding vector or high-resolution images can be viewed in the result files. The result files will be released together with the original data.
References
- Li H,Durbin R.Fast and accurate short read alignment with Burrows-Wheeler Transform.Bioinformatics, 2009 25:1754-60
- Danecek, Petr, et al. Twelve years of SAMtools and BCFtools. Gigascience 10.2 (2021): giab008.
- McKenna A, Hanna M, Banks E, Sivachenko A, et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010 20:1297-303
- Cingolani P, Platts A, Wang le L, et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3., Fly (Austin). 2012 Apr-Jun;6(2):80-92.
- Manta: rapid detection of structural variants and indels for germline and cancer sequencing applications. Bioinformatics, 32, 1220-1222.
- Boeva V, Popova T, Bleakley K, Chiche P, Cappo J, Schleiermacher G, Janoueix-Lerosey I, Delattre O, Barillot E. (2012) Control-FREEC: a tool for assessing copy number and allelic content using next generation sequencing data. Bioinformatics. Bioinformatics, 2012, 28(3):423-5. PubMed PMID: 22155870.
- Sensitive protein alignments at tree-of-life scale using DIAMOND
- Yangyang DENG, Jianqi LI, Songfeng WU, et al. Integrated nr Database in Protein Annotation System and Its Localization. Computer Engineering, 2006 32(5):71_74
- Yangyang DENG, Jianqi LI, Songfeng WU, et al. Integrated nr Database in Protein Annotation System and Its Localization. Computer Engineering, 2006 32(5):71_74
- Michael Ashburner, Catherine A. Ball, Judith A. Blake, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet, (25): 25_29
- Roman L.Tatusov, Michael Y.Galperin, Darren A.Natale, et al. The COG database: a tool for genome_scale analysis of protein functions and evolution. Nucleic Acids Res, 2000_7_1 28(1):33_6
- Minoru Kanehisa, Susumu Goto, Shuichi Kawashima, et al. The KEGG resource for deciphering the genome. Nucleic Acids Res 2004, (32):D277_D280