前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >高通量测序如何寻找T-DNA插入的位置

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

作者头像
生信技能树
发布2018-03-29 16:25:09
17.4K0
发布2018-03-29 16:25:09
举报
文章被收录于专栏:生信技能树

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

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

用efecth下载参考基因组

代码语言:javascript
复制
mkdir -p refsefetch -db=nuccore -format=fasta -id=AF086833 | seqret --filter -sid AF086833 > refs/AF086833.fa

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

代码语言:javascript
复制
# 提取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安装

代码语言:javascript
复制
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为前缀的几个文件。

第三步:回贴到参考序列。所用工具为 bwasamtools

代码语言:javascript
复制
# 建立索引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统计比对的总体情况:

代码语言:javascript
复制
$ 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%序列是什么情况?我们可以通过管道的方式进行简单的统计

代码语言:javascript
复制
$ 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

其中大部分序列是 77141,也就是说两条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。

代码语言:javascript
复制
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安装。

代码语言:javascript
复制
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的可能位置。

代码语言:javascript
复制
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

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

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2018-03-11,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信技能树 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
相关产品与服务
云直播
云直播(Cloud Streaming Services,CSS)为您提供极速、稳定、专业的云端直播处理服务,根据业务的不同直播场景需求,云直播提供了标准直播、快直播、云导播台三种服务,分别针对大规模实时观看、超低延时直播、便捷云端导播的场景,配合腾讯云视立方·直播 SDK,为您提供一站式的音视频直播解决方案。
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档