高通量测序如何寻找T-DNA插入的位置

为了解基因组存在T-DNA插入时，即基因组构成为AC而样本基因组为ABC的情况得到的测序结果在序列比对的时候的可能情况，因此需要先要使用模拟数据进行探索。

第一步：构建参考序列和实际序列。这一部分会用到 samtools, emboss和 entrez-direct, 都可以通过conda安装

用efecth下载参考基因组

  1. mkdir -p refs

  2. efetch -db=nuccore -format=fasta -id=AF086833 | seqret --filter -sid AF086833 > refs/AF086833.

  3. fa

从参考基因组提取其中部分序列用作参考序列，而下载的参考基因组则被当成实际的基因组

  1. # 提取1~5000, 8000~

  2. cat refs/AF086833.fa | seqret -filter -sbegin 1 -send 5000 > part1.fa

  3. cat refs/AF086833.fa | seqret -filter -sbegin 8000 > part2.fa

  4. # 合并

  5. cat part1.fa part2.fa| union -filter > refs/ref.fa

第二步：模拟测序结果。这一步用到 dwgsim，也可以用 conda安装

  1. mkdir data

  2. dwgsim -e 0.02 -E 0.02 -d 350 -1 100 -2 100 -s 50 -r 0 -R 0 -N 10000 -c 0 refs/AF086833.fa data/data

解释dwgsim的参数， -e和 -E为测序仪的系统错误率, -d表示文库大小, -1和 -2表示短读长度(这里就是文库大小350bp，PE100), 而 -s则表示文库大小的波动情况， -r和 -R表示基因组的突变率, -N表示输出的短读数， -c表示输出数据类型(0为illumina, 1为SOLiD,2为Ion Torrent)。最后会在data文件下生成以data为前缀的几个文件。

第三步：回贴到参考序列。所用工具为 bwa和 samtools

  1. # 建立索引

  2. bwa index refs/ref.fa

  3. bwa mem refs/ref.fa

  4. # 比对

  5. mkdir align

  6. bwa mem refs/ref.fa data/data.bwa.read1.fastq.gz data/data.bwa.read2.fastq.gz| samtools sort > align/data.bwa.bam

  7. samtools index align/data.bwa.bam

第四步：使用IGV和samtools探索比对结果. samtools是处理SAM/BAM格式的常用工具，而IGV则是可视化利器。首先用samtools的flagstat统计比对的总体情况：

  1. $ samtools flagstat align/data.bwa.bam

  2. 20000 + 0 in total (QC-passed reads + QC-failed reads)

  3. 0 + 0 secondary

  4. 0 + 0 supplementary

  5. 0 + 0 duplicates

  6. 15906 + 0 mapped (79.53% : N/A)

  7. 20000 + 0 paired in sequencing

  8. 10000 + 0 read1

  9. 10000 + 0 read2

  10. 15662 + 0 properly paired (78.31% : N/A)

  11. 15662 + 0 with itself and mate mapped

  12. 244 + 0 singletons (1.22% : N/A)

  13. 0 + 0 with mate mapped to a different chr

  14. 0 + 0 with mate mapped to a different chr (mapQ>=5)

不难大部分序列(~80%)都是正确成对(properly paired)，其中properly paired的解释为"0x2 PROPER_PAIR .. each segment properly aligned according to the aligner"，也就是两个序列都能在基因组上找到自己的位置，最常见的两类flags就是"83,163"和"99和147",也就是和参考序列反向互补。

flags为83和163的结果

那么余下的20%序列是什么情况？我们可以通过管道的方式进行简单的统计

  1. $ samtools view -F 0x2 align/data.bwa.bam | cut -f 2 | sort | uniq -c | sort -k1,1nr

  2. 1925 141

  3. 1925 77

  4. 65 137

  5. 65 69

  6. 62 121

  7. 62 181

  8. 59 117

  9. 59 185

  10. 58 133

  11. 58 73

其中大部分序列是 77和 141，也就是说两条reads都没有比对到参考基因组上， 也就是SAM格式中的第3，6，7列为"*",第4,5,8,9列表示为"0"

141和77表示完全没有比对

剩下的"69,137","117,185"和"73,133","121,181"表示两条reads中只有一条,即flags为137,185和73,121的reads能比对到参考基因组上。

关于flags的含义，可以使用网页版 https://broadinstitute.github.io/picard/explain-flags.html 查询

其中一条read比对到参考序列

如果统计这些单边比对reads的位置信息，就会发现他们的位置是在4651~5214, 也就缩小搜索区间，因为通过IGV你会发现区间刚好存在一个breakpoint，所有双端联配在这里都出现不同程度的soft-clip。

  1. samtools view -b -F 0x2 align/data.bwa.bam | samtools view -b -G 141 | samtools view -G 77 | cut -f 4 | sort | head -n2

  2. samtools view -b -F 0x2 align/data.bwa.bam | samtools view -b -G 141 | samtools view -G 77 | cut -f 4 | sort | tail -n2

5000bp处就是插入位置

第五步：组装20%的不完美比对序列。这一步使用velvet组装工具,因为用起来比较容易，而且可以用bioconda安装。

  1. samtools view -b -F 0x2 align/data.bwa.bam | samtools sort -n | samtools fastq -1 read_1.fq -2 read_2.fq -

  2. velveth velvet31 31 -fastq -separate -shortPaired read_1.fq read_2.fq

  3. velvetg velvet31 -exp_cov auto -ins_length 150

  4. # -ins_length表示两个reads的间隔平均距离

最后会在velvet31文件夹下生成 contigs.fa,这里面的N50肯定是看不了的，我们只是需要一个比较长一点的序列而已。

第六步：使用BLAST找到可能的位点。建立索引数据库，然后搜索组装的 contigs.fa的可能位置。

  1. cd refs

  2. # 建立索引

  3. makeblastdb -dbtype nucl -in ref.fa

  4. # 搜索

  5. blastn -query ../velvet31/contigs.fa -db ref.fa -outfmt 8

BLASTN结果

以上仅仅使用了模拟数据的方式验证了方案的可行性，实际情况会比较复杂

参考文献

- "Genetic characterization of T-DNA insertions in the genome of the Arabidopsis thaliana sumo1/2 knock-down line

- Illumina Sequencing Technology as a Method of Identifying T-DNA Insertion Loci in ActivationTagged Arabidopsis thaliana Plants

- 【直播】我的基因组（十五）:提取未比对的测序数据

赞 (0)

相关推荐