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

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

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

用efecth下载参考基因组

mkdir -p refsefetch -db=nuccore -format=fasta -id=AF086833 | seqret --filter -sid AF086833 > refs/AF086833.fa

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

# 提取1~5000, 8000~cat refs/AF086833.fa | seqret -filter -sbegin 1 -send 5000 > part1.facat refs/AF086833.fa | seqret -filter -sbegin 8000 > part2.fa# 合并cat part1.fa part2.fa| union -filter > refs/ref.fa

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

mkdir datadwgsim -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

# 建立索引bwa index refs/ref.fabwa mem refs/ref.fa# 比对mkdir alignbwa mem refs/ref.fa data/data.bwa.read1.fastq.gz data/data.bwa.read2.fastq.gz| samtools sort > align/data.bwa.bamsamtools index align/data.bwa.bam

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

$ samtools flagstat align/data.bwa.bam20000 + 0 in total (QC-passed reads + QC-failed reads)0 + 0 secondary0 + 0 supplementary0 + 0 duplicates15906 + 0 mapped (79.53% : N/A)20000 + 0 paired in sequencing10000 + 0 read110000 + 0 read215662 + 0 properly paired (78.31% : N/A)15662 + 0 with itself and mate mapped244 + 0 singletons (1.22% : N/A)0 + 0 with mate mapped to a different chr0 + 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%序列是什么情况？我们可以通过管道的方式进行简单的统计

$ samtools view -F 0x2 align/data.bwa.bam | cut -f 2 | sort | uniq -c | sort -k1,1nr1925 1411925 7765 13765 6962 12162 18159 11759 18558 13358 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。

samtools view -b -F 0x2 align/data.bwa.bam | samtools view -b -G 141 | samtools view -G 77 | cut -f 4 | sort | head -n2samtools 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安装。

samtools view -b -F 0x2 align/data.bwa.bam | samtools sort -n | samtools fastq -1 read_1.fq -2 read_2.fq -velveth velvet31 31 -fastq -separate -shortPaired read_1.fq read_2.fqvelvetg velvet31 -exp_cov auto -ins_length 150# -ins_length表示两个reads的间隔平均距离

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

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

cd refs# 建立索引makeblastdb -dbtype nucl -in ref.fa# 搜索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

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