Screening method of variation sites in fragmented DNA re-sequencing data

Abstract

The invention discloses a method for screening variation sites in fragmented DNA re-sequencing data, and belongs to the technical field of biology. In order to solve the problem that a method for efficiently and accurately screening variation sites in fragmented DNA re-sequencing data lacks in the prior art, the invention provides a method for screening variation sites in fragmented DNA re-sequencing data, according to the screening method, when-fq is converted into a. Sai file in sequence alignment,-lmin-size-n 0.04 is selected as a parameter; when-. Sai is converted into. Sam in sequence alignment,-o = 1000 is selected as a parameter; when variation sites are screened, the threshold value of the selection parameter MQ is 36. Compared with a traditional variation site screening method, the screening method provided by the invention has the advantages that the efficient comparison rate is ensured, the operation time of differential site screening is remarkably shortened, the accuracy of less selection results is improved, and the method can be applied to screening of variation sites in fragmented DNA re-sequencing data.

Description

Technical Field

[0001] The present invention belongs to the field of biotechnology, and in particular relates to a method for screening variant sites in fragmented DNA resequencing data. Background Art

[0002] Genomic DNA records all the genetic information of living organisms. This genetic information can promote research in many aspects such as species identification, genetic diversity assessment, disease monitoring and ecological adaptation, and provide a strong theoretical basis for scientists to better understand and protect wild animals and create a good ecological environment. Materials that can only obtain genomic DNA from some historical specimens or keratinized materials do not have a complete cell structure, so the nuclear genetic material inside them will be attacked by various enzymes, causing them to continuously degrade. Studies have shown that genetic material can be preserved for tens of thousands of years in specific environments such as hair, bones, teeth, and soft tissues. Through these genetic materials, we can understand the genetic diversity of a certain group, reveal the changes in the ecosystem of a certain area, and even the whole process of some animal groups from prosperity to decline and eventually to extinction in history. At present, the origin research of various known species is also closely related to fragmented DNA.

[0003] The emergence of Sanger sequencing technology (first-generation sequencing technology) marks the visualization of genetic information as a reality, but there are technical problems such as throughput and read length limitations. The emergence and continuous improvement of second-generation sequencing technology (NGS) means that large-scale sequencing of small DNA fragments has become possible. However, NGS can generate billions of DNA sequences, which requires accurate alignment of these fragments with the corresponding reference genome, and it takes a lot of time to find the different locations (i.e., variant sites).

[0004] Therefore, those skilled in the art are eager to develop an efficient and accurate method for screening variant sites in fragmented DNA, especially small fragmented DNA resequencing data. Summary of the invention

[0005] The present invention solves the problem that there is a lack of efficient and accurate method for screening variant sites in fragmented DNA resequencing data in the prior art, and provides a method for screening variant sites in fragmented DNA resequencing data.

[0006] One of the purposes of the present invention is to provide a method for screening variant sites in fragmented DNA resequencing data, the screening method comprising the following steps: aligning the fragmented DNA resequencing data to be tested to a reference genome sequence, performing sequence alignment using BWA software, filtering variant sites on the file results obtained by the above sequence alignment, and screening to obtain variant sites on the fragmented DNA genome to be tested.

[0007] In a preferred embodiment of the present invention, the sequence alignment steps are as follows: in the first step of the sequence alignment, the fq file is first converted into a .sai file, and the command used is bwa aln-l min-size-n 0.04-t 5$.ref.fa$.1.fq.gz-f$.1.sai&bwa aln-l min-size-k 2-t 5$.ref.fa$.2.fq.gz-f$.2.sai, and the effect is evaluated by changing the parameters Rm, Rm0 and Tc; in the second step of the sequence alignment, the above-mentioned .sai file is converted into a .sam file, and the command used is bwa sampe-o 1000$.ref.fa$.1.sai$.2.sai$.1.fq.gz$.2.fq.gz>out.sam, and the effect is evaluated by changing the parameters Rm, Rm0 and Tc.

[0008] In a preferred embodiment of the present invention, in the first step of the sequence alignment, the parameter -l is the seed length of the alignment, min-size is the minimum read length, -n is the error rate, -t is the number of threads, $.ref.fa is the reference genome, $.1.fq.gz / $.2.fq.gz is the input fq file 1 / 2, and $.1.sai / $.2.sai is the corresponding output file 1 / 2.

[0009] In a preferred embodiment of the present invention, the parameter -o in the second step of the sequence alignment is the number of times the alignment is repeated.

[0010] In a preferred embodiment of the present invention, the parameter Rm is the alignment rate, specifically refers to the ratio of the reads generated by sequencing that are aligned to the reference genome to the total reads.

[0011] In a preferred embodiment of the present invention, the parameter Rm0 is the multiple alignment rate, specifically refers to the ratio of the number of reads of the same read aligned to different reference sequences to the number of reads on the pair.

[0012] In a preferred embodiment of the present invention, the parameter Tc is the running time, and c is the number of running threads.

[0013] In a preferred embodiment of the present invention, the command used in the variant site filtering step is gatkVariantFilter-R$.ref.fa-V$.snp.vcf--filter-expression"MQ<36.0||QD<3.0||FS>60.0||SOR>3.0||MQRankSum<-12.5||ReadPosRankSum<-8.0"--filter-name"Filter"--cluster-window-size 10--cluster-size 3--missing-values-evaluate-as-failing-O MQ36.$.snp.vcf.

[0014] In a preferred embodiment of the present invention, in the variant site filtering step, parameter -V is the input file, --filter-expression is the filtering parameter, --cluster-window-size is the window length, --cluster-size is the number of sites allowed in each window, and -O is the output file.

[0015] Beneficial effects of the present invention: The present invention provides a method for screening variant sites in fragmented DNA resequencing data. When -fq is converted to .sai file in sequence alignment, the minimum reading length is used as the alignment seed length (-l), and the threshold of the error rate of the alignment (-n) is changed from the traditional "1024-n 0.03" parameter to "20-n 0.04"; when -.sai is converted to .sam file in sequence alignment, the number of repeated alignments of a segment of reading (-o) is changed from 100000 to 1000; when screening variant sites, according to the characteristics of the fragmented DNA resequencing data reads, the threshold of the parameter MQ is changed from 40 to 36.

[0016] The method for screening variant sites in fragmented DNA resequencing data provided by the present invention, compared with the traditional method for screening variant sites, significantly shortens the running time of the comparison while ensuring an efficient comparison rate, and can maximize the screening efficiency of variant sites in fragmented DNA resequencing data; compared with the traditional screening method, the method for screening variant sites in fragmented DNA resequencing data provided by the present invention is more suitable for screening variant sites in fragmented DNA, especially small fragmented DNA resequencing data, significantly improves the screening efficiency and accuracy of differential sites, significantly increases the number of variant sites obtained by screening, and provides more basic sites for subsequent analysis; and the screening method provided by the present invention has little impact on high-quality DNA sequencing data, has a wide range of applications, and can be applied to the screening of variant sites in fragmented DNA resequencing data. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 It is the comparison efficiency diagram of S1-S5 samples under different parameters; A is the comparison rate (Rm) result diagram, B is the multiple comparison rate (Rm0) result diagram, and C is the running time (Tc) result diagram;

[0018] Figure 2 The results of the effect of different -o values ​​on the comparison efficiency are shown in Figure 1. A is the multiple comparison rate (Rm0) result diagram, and B is the running time (Tc) result diagram.

[0019] Figure 3 A is a characteristic comparison diagram of SNP sites in hair and skin tissue samples; A is the MQ distribution diagram of SNP sites, and B is the distribution diagram of reading length when MQ is 37 or 60. DETAILED DESCRIPTION

[0020] Those skilled in the art can refer to the content of this article and appropriately improve the process parameters. It should be particularly noted that all similar replacements and modifications are obvious to those skilled in the art, and they are all considered to be included in the present invention. The method and application of the present invention have been described through preferred embodiments, and relevant personnel can obviously modify or appropriately change and combine the methods and applications described herein without departing from the content and scope of the present invention to implement and apply the technology of the present invention.

[0021] In order to make the purpose, technical scheme and advantages of the present invention more clear, the present invention is further described in detail below in conjunction with specific embodiments and the accompanying drawings of the specification. The test methods used in the following examples are conventional methods unless otherwise specified, and the materials, reagents, methods and instruments used are conventional materials, reagents, methods and instruments in the art unless otherwise specified, and can be obtained by those skilled in the art through commercial channels.

[0022] Example 1: Alignment of sequencing reads of fragmented DNA to a reference genome

[0023] The research objects of this embodiment include tortoise shell (S1), horn (S2), pangolin scale (S3), leopard cat (S4) and arctic fox hair (S5), and 100000000 readings are extracted from each sample as the sample to be tested. The data comes from the NCBI database, the project number is PRJNA1103004, and the data numbers corresponding to samples S1-S5 are SRR28762820, SRR28762814, SRR28762863, SRR28762865, SRR28762854, and the corresponding reference genome numbers are GCF_016161935.1, GCF_019923935.1, GCF_014570535.1, GCF_016509475.1 and GCF_018345385.1. Align the DNA resequencing data of the above-mentioned sample to be tested to the reference genome sequence, and use BWA software for sequence alignment. The specific command for -fq to convert to .sai file is as follows:

[0024] bwa aln-l 20-n 0.04-t 5$.ref.fasta$.1.fq.gz-f$.1.sai;

[0025] bwa aln-l 20-n 0.04-t 5$.ref.fasta$.2.fq.gz-f$.2.sai;

[0026] The above .sai file is converted into a .sam file by command "bwa sampe-o 1000-r'@RG\tID:$ID\tSM:$1\tPL:UNKNOWN'$.ref.fa$.1.sai$.2.sai$.1.fq.gz$.2.fq.gz>$.out.sam", and Rm, Rm0 and Tc are counted; next, the .sai file generated by "-l 20-n 0.04" continues to discuss the value range of -o, and -o is 1000 / 10000 / 100000 / 1000000 respectively. The specific commands are as follows:

[0027] bwa sampe-o 1000-r'@RG\tID:$ID\tSM:$1\tPL:UNKNOWN'$.ref.fa$.1.sai$.2.sai$.1.fq.gz$.2.fq.gz>$.out.sam;

[0028] bwa sampe-o 10000-r'@RG\tID:$ID\tSM:$1\tPL:UNKNOWN'$.ref.fa$.1.sai$.2.sai$.1.fq.gz$.2.fq.gz>$.out.sam;

[0029] bwa sampe-o 100000-r'@RG\tID:$ID\tSM:$1\tPL:UNKNOWN'$.ref.fa$.1.sai$.2.sai$.1.fq.gz$.2.fq.gz>$.out.sam;

[0030] bwa sampe-o 1000000-r'@RG\tID:$ID\tSM:$1\tPL:UNKNOWN'$.ref.fa$.1.sai$.2.sai$.1.fq.gz$.2.fq.gz>$.out.sam.

[0031] Example 2: Characteristics and screening of variant sites

[0032] In this example, the hair of an Arctic fox is selected as the sample to be tested (B1-H), and its skin tissue is used as the control group (B1-S). The data of the sample to be tested comes from the NCBI database, project number: PRJNA1103004, and the data numbers of samples B1-H and BI-S are SRR28762875 and SRR28762874, respectively. About 75G of data are randomly selected from the above B1-H and BI-S samples for SNP calling, and the reference genome number is GCF_018345385.1. The specific command is as follows:

[0033] bwa aln-l 20-n 0.04-t 5ref.fasta$.1.fq.gz-f$.1.sai;

[0034] bwa aln-l 20-n 0.04-t 5ref.fasta$.2.fq.gz-f$.2.sai;

[0035] bwa sampe-o 1000-r'@RG\tID:$ID\tSM:$1\tPL:UNKNOWN'ref.fa$.1.sai$.2.sai$.1.fq.gz$.2.fq.gz>$.out.sam;

[0036] samtools view-Shb$.out.sam>$.out.bam;

[0037] sambamba view-t 5-hf bam. / $.out.bam-o$.filtered_output.bam;

[0038] sambamba sort -t 5 --tmpdir=$PWD -p. / $.filtered_output.bam -o $.out.sort.bam;

[0039] sambamba markdup -t 5 --tmpdir=$PWD -p $.out.sort.bam o $.ut.sort.markdup.bam;

[0040] sambamba sort -t 5 --tmpdir=$PWD -p. / $.out.sort.markdup.bam -o $.out.sort.markdup1.bam;

[0041] gatk HaplotypeCaller --native-pair-hmm-threads 5 -R ref.fa -I $.out.sort.markdup1.bam -O $.raw.vcf;

[0042] gatk SelectVariants -R ref.fna -V. / $.raw.vcf --select-type-to-include SNP -O $.snp.raw.vcf;

[0043] gatk VariantFiltration -R ref.fna -V. / $.snp.raw.vcf --filter-expression "MQ<40.0||QD<2.0||FS>60.0||SOR>3.0||MQRankSum<-12.5||ReadPosRankSum<-8.0" --filter-name "Filter" --cluster-window-size 10 --cluster-size 3 --missing-values-evaluate-as-failing -O $.snp.filter.vcf.tmp;

[0044] gatk SelectVariants -R ref.fna -V. / $.snp.filter.vcf.tmp --exclude-filtered -O $.MQ40_snp.filter.vcf。

[0045] Comparative Example 1:

[0046] The difference between this comparative example and Example 1 is that the specific command of -fq to convert .sai file is as follows:

[0047] bwa aln-l 1024-n 0.03-t 5$.ref.fasta$.1.fq.gz-f$.1.sai;

[0048] bwa aln-l 1024-n 0.03-t 5$.ref.fasta$.2.fq.gz-f$.2.sai.

[0049] The other steps and parameters are the same as those in Example 1.

[0050] Comparative Example 2:

[0051] The difference between this comparative example and Example 2 is that the threshold value of the parameter MQ is 37. The other steps and parameters are the same as those in Example 1.

[0052] Effect experiment:

[0053] (1) The present invention performs statistics of Rm, Rm0 and Tc on the result files after -fq is converted to .sai files in Example 1 and Comparative Example 1, and compares and analyzes the processing results of different comparison parameters. The results are as follows Figure 1 As shown, the comparison efficiency (Rm) of the five samples to be tested (S1-S5) under the two comparison strategies of Example 1 and Comparative Example 1 is almost the same (p=0.9983), and the Rm values ​​are 66.42%±21.51% and 65.25%±22.12% respectively. Figure 1 A); There was no significant difference in the multiple comparison rates between Example 1 and Comparative Example 1 under the two comparison strategies (p = 0.9993) ( Figure 1 B); There is a significant difference in the running time under the two comparison strategies of Example 1 and Comparative Example 1. The running time required for Comparative Example 1 is 3.19±0.80 times that required for Example 1 (p=0.0000) ( Figure 1 C). It can be seen that the method for screening variant sites in fragmented DNA resequencing data provided by the present invention significantly shortens the running time of the alignment while ensuring a high alignment rate compared to the traditional screening method.

[0054] (2) The present invention uses the solution 20_0.04 described in Example 1 to convert -fq to .sai file as a premise, and analyzes the value of -o. The results are as follows Figure 2 As shown, different -o values ​​have no effect on the comparison rate ( Figure 2A), the alignment rates were all 65.25% ± 22.12%; in terms of multiple alignment rate (Rm0), when the -o value increased from 1000 to 10000, 100000 and 1000000, the Rm0 was only reduced by 0.60% ± 0.48%, 1.42% ± 1.02% and 1.45% ± 1.04% on average, respectively, and the differences between the groups were significant (P1000vs10000=0.0514, P1000vs100000=0.0365, P1000vs1000000=0.0358)( Figure 2 A). Compared with the experimental group with -o value of 1000, the running time (Tc) of the experimental groups with -o value of 10000, 100000, and 100000 increased by 4.43±1.65 (P1000vs10000=0.0000), 24.87±21.36 (P1000vs10000=0.0000), and 34.49±29.45 (P1000vs10000=0.0000) ( Figure 2 B). It can be seen that when the -o value is 1000, the required running time is the shortest, which can maximize the screening efficiency of variant sites in fragmented DNA resequencing data.

[0055] (3) The present invention further explored the characteristics of SNP typing sites in hair (B1-H) and skin (B1-S) sample DNA. Figure 3 As shown, 5145653 and 5285155 raw variant sites (Nraw) were obtained from the data of B1-H and B1-S, respectively, and the difference in the number of variant sites was not significant. However, in the results obtained using the conventional commands described in Comparative Example 1, B1-S retained 3497794 variant sites (N40), while B1-H only retained 1534471 variant sites.

[0056] The present invention analyzes the SNP sites of hair (B1-H) and skin (B1-S) samples, and the results are as follows Figure 3 As shown in Figure 2, the MQ values ​​of the variant sites in the hair sample (B1-H) were concentrated at 37, while the MQ values ​​of the skin samples were concentrated at 60 ( Figure 3 A). It can be seen that when the threshold of the parameter MQ of the conventional filtering in Comparative Example 2 is 37, a large number of sites of the hair sample (B1-H) will be filtered out, thereby causing inaccurate screening results.

[0057] Therefore, after adjusting the threshold of the conventional filtering parameter MQ to 36, the number of SNP sites in the B1-H sample increased by 166.48% to 2554541; and it was found that as long as the length of the read was less than the maximum read, the MQ quality score might be deducted accordingly, which eventually led to the MQ values ​​of these perfectly aligned reads being concentrated around 37 ( Figure 3 B) Therefore, for these highly fragmented DNA reads, setting the MQ threshold below 37 can significantly increase the number of sites when calling SNPs.

[0058] It can be seen that the screening method for variant sites in fragmented DNA resequencing data provided by the present invention is more suitable for fragmented DNA, especially the screening of variant sites in small-fragment DNA resequencing data, compared with traditional screening methods, and significantly shortens the running time of the comparison, has higher screening efficiency, and the screening results are more accurate.

[0059] The contents not described in detail in the specification of the present invention are well-known technologies for those skilled in the art. Although the present invention has been disclosed as above with preferred embodiments, it is not intended to limit the present invention. Anyone familiar with this technology can make various changes and modifications without departing from the spirit and scope of the present invention. Therefore, the protection scope of the present invention should be based on the definition of the claims.

Claims

1. A method for screening variant sites in fragmented DNA resequencing data, characterized in that: The screening method comprises the following steps: aligning the fragmented DNA resequencing data to be tested to the reference genome sequence, performing sequence alignment using BWA software, filtering the variant sites on the file results obtained by the sequence alignment, and screening to obtain the variant sites on the fragmented DNA genome to be tested.

2. The screening method according to claim 1, characterized in that The sequence alignment steps are as follows: in the first step of the sequence alignment, the fq file is first converted into a .sai file, and the command used is bwa aln-l min-size-n 0.04-t 5$.ref.fa$.1.fq.gz-f$.1.sai&bwa aln-l min-size-k 2-t 5$.ref.fa$.2.fq.gz-f$.2.sai, and the effect is evaluated by changing the parameters Rm, Rm0 and Tc; in the second step of the sequence alignment, the above .sai file is converted into a .sam file, and the command used is bwa sampe-o 1000$.ref.fa$.1.sai$.2.sai$.1.fq.gz$.2.fq.gz>out.sam, and the effect is evaluated by changing the parameters Rm, Rm0 and Tc.

3. The screening method according to claim 2, characterized in that In the first step of the sequence alignment, the parameter -l is the seed length of the alignment, min-size is the minimum read length, -n is the error rate, -t is the number of threads, $.ref.fa is the reference genome, $.1.fq.gz / $.2.fq.gz is the input fq file 1 / 2, and $.1.sai / $.2.sai is the corresponding output file 1 / 2.

4. The screening method according to claim 2, characterized in that In the second step of the sequence alignment, the parameter -o is the number of times the alignment is repeated.

5. The screening method according to claim 2, characterized in that The parameter Rm is the alignment rate, specifically refers to the ratio of the reads generated by sequencing that are aligned to the reference genome to the total reads.

6. The screening method according to claim 2, characterized in that The parameter Rm0 is the multiple alignment rate, specifically refers to the ratio of the number of reads of the same read aligned to different reference sequences to the number of reads in the pair.

7. The screening method according to claim 2, characterized in that The parameter Tc is the running time, and c is the number of running threads.

8. The screening method according to claim 1, characterized in that The command used in the variant site filtering step is gatk VariantFiltration-R$.ref.fa-V$.snp.vcf--filter-expression"MQ<36.0||QD<3.0||FS>60.0||SOR>3.0||MQRankSum<-12.5||ReadPosRankSum<-8.0"--filter-name"Filter"--cluster-window-size 10--cluster-size 3--missing-values-evaluate-as-failing-O MQ36.$.snp.vcf.

9. The screening method according to claim 8, characterized in that In the variant site filtering step, the parameter -V is the input file, --filter-expression is the filtering parameter, --cluster-window-size is the window length, --cluster-size is the number of sites allowed in each window, and -O is the output file.

Images