Loading [MathJax]/jax/output/CommonHTML/config.js
前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
圈层
工具
发布
首页
学习
活动
专区
圈层
工具
MCP广场
社区首页 >专栏 >高通量测序如何寻找T-DNA插入的位置

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

作者头像
生信技能树
发布于 2018-03-29 08:25:09
发布于 2018-03-29 08:25:09
17.5K00
代码可运行
举报
文章被收录于专栏:生信技能树生信技能树
运行总次数:0
代码可运行

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

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

用efecth下载参考基因组

代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
mkdir -p refsefetch -db=nuccore -format=fasta -id=AF086833 | seqret --filter -sid AF086833 > refs/AF086833.fa

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

代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
# 提取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
代码运行次数:0
运行
AI代码解释
复制
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
代码运行次数:0
运行
AI代码解释
复制
# 建立索引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
代码运行次数:0
运行
AI代码解释
复制
$ 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
代码运行次数:0
运行
AI代码解释
复制
$ 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
代码运行次数:0
运行
AI代码解释
复制
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
代码运行次数:0
运行
AI代码解释
复制
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
代码运行次数:0
运行
AI代码解释
复制
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 删除。

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

评论
登录后参与评论
暂无评论
推荐阅读
编辑精选文章
换一批
Greenplum 7 新特性整理
参考:https://www.xmmup.com/zaidockerzhongkuaisutiyangreenplum-7-0-0.html
AiDBA宝典
2023/10/16
1.4K0
Greenplum 7 新特性整理
Postgresql10分区表range实例
创建规则 https://www.postgresql.org/docs/10/sql-createtable.html
mingjie
2022/05/12
3680
Postgresql分区表大量实例与分区建议(LIST / RANGE / HASH / 多级混合分区)
5.11.6. Best Practices for Declarative Partitioning
mingjie
2022/09/26
7.1K0
【动手实践】Oracle 12.2 新特性:自动的列表分区创建
2017年来了,我们要启动新的学习征程了。在过去我们一直思考,什么样的内容能够更帮助大家了解和学习到有用的知识? 这个『动手实践』栏目就是这样一个改进和尝试吧,一个小小的范例,几分钟的线上实践(感谢Oracle),就能帮助大家熟悉一个知识点,几个重要的命令。如此是否会有不一样的体验?试一试吧。 ---- 在Oracle Database 12.2 之前,如果使用列表分区,当插入的数据超过了分区列表值设定,则会抛出异常;而如果存在大量的列表值需要定义,则可能需要一一设置。 在12.2引入的新特性中 - Au
数据和云
2018/03/06
1.2K0
【动手实践】Oracle 12.2 新特性:自动的列表分区创建
zabbix表分区(适用于zabbix2.0.x,zabbix2.2.x和zabbix2.4.x)[推荐]
本文主要介绍了zabbix进行数据库表分区的方法: 在系统监控中,zabbix已经代替了nagios+cacti,zabbix以其良好的图形展示和高度自定义赢得了很多运维人员的喜爱。但是由于在工作中,zabbix跑的时间过长(我们公司跑了将近3年),web页面经常卡顿,监控数据有时很难插入数据库,且数据库队列经常性卡死,经过查看,发现mysql的数据量高达83G,急需瘦身,于是有了此文。 步骤: 修改表结构: use zabbix; Alter table history_text drop
小小科
2018/05/04
7920
【Z投稿】大规模数据库监控的Zabbix玩法详谈
10多年MySQL大规模数据库运维经验+8年Zabbix使用经验。本次峰会演讲和workshop《大规模数据库监控的Zabbix玩法》,讲述海量数据库实例的监控,介绍zabbix的安装、部署、优化,以及数据库自动化运维。
Zabbix
2021/02/03
7230
Database of Databases 搜索引擎的妙用
最近在用 Database of Databases 去查询TiDB 、PolarDB、SequoiaDB、OceanBase的相关资料并做了一些简单的对比。比如像下面那个表格一样,可以对这四个数据库的基本信息进行对比。
哒呵呵
2021/07/19
7920
mysql分区语句
要是分区数比现有的分区数多的话,只能使用 ADD来添加分区数.下面就表示增加了6个分区数
全栈程序员站长
2022/08/11
12.7K0
Zabbix MySQL MariaDB 数据库分表
Zabbix 数据库在没有使用分区分表功能,默认使用Housekeeping(管家功能)进行删除历史数据和趋势历史记录,如果zabbix数据库使用了分区分表功能需要把Housekeeping(管理功能)关闭。Housekeeping功能监控数据量少可以使用,但监控数据量多每次执行删除旧数据会降低MySQL数据库性能,并且还会产生很多空间碎片。经常会出现警报" Zabbix housekeeper processes more than 75% busy"的告警。(zabbix_server.conf配置文件两个参数进行历史记录数据删除:间隔多久删除一次,默认单位小时HousekeepingFrequency=1,一次删除多少数据,默认单位行MaxHousekeeperDelete=5000)。
Kevin song
2021/09/15
2.1K0
Hive 各版本关键新特性(Key New Feature)介绍
开源世界里的代码受社区推动和极客文化的影响,变化一直都很快。这点在 hadoop 生态圈里表现尤为突出,不过这也与 hadoop 得到业界的广泛应用以及各种需求推动密不可分(近几年大数据、云计算被炒烂
用户1177713
2018/02/24
2.6K0
Hive 各版本关键新特性(Key New Feature)介绍
SQL Server通过整理索引碎片和重建索引提高速度
本文章转载:http://database.51cto.com/art/201108/282408.htm
跟着阿笨一起玩NET
2018/09/19
4.5K0
Zabbix 分区优化
文章主要介绍了如何通过修改MySQL配置文件中的PARTITION参数来达到提升性能的目的。具体来说,可以通过设置合适的PARTITION值来减少分区数,从而提高查询速度。同时,建议使用InnoDB引擎,并合理配置参数,以达到最佳性能。
腾讯云TStack
2017/11/09
3.7K0
ubuntu分区方案(合并分区)
一句话概括:Ubuntu系统在一个硬盘上只支持最多4个 Primary 分区或3个 Primary 分区加1个 Extended 分区。Extended 分区下面可以有多个 Logical 分区。
全栈程序员站长
2022/07/29
5K0
Hive & Performance 学习笔记
注:本文来源于 Hortonworks 的 Adam Muise 在 July 23 2013 日的 Toronto Hadoop User Group 大会上的一次演讲, 本文只是稍作增删、整理
用户1177713
2018/02/24
1.5K0
Hive & Performance 学习笔记
Apache Spark 2.2.0 中文文档 - Spark SQL, DataFrames and Datasets Guide | ApacheCN
本文介绍了基于Spark的SQL编程的常用概念和技术。首先介绍了Spark的基本概念和架构,然后详细讲解了Spark的数据类型和SQL函数,最后列举了一些Spark在实际应用中的例子。
片刻
2018/01/05
26.5K0
Hive SQL 语法大全,宇宙最强整理,建议收藏
LOCATION 是指定外部表的存储路径,MANAGEDLOCATION 是指定管理表的存储路径(hive 4.0.0 才支持),官方建议默认就行,让所有的表都在一个根目录下。
kk大数据
2020/11/03
7.3K0
SQL server 2005 切换分区表
如转载,请注明出处:http://blog.csdn.net/robinson_0612/archive/2009/11/10/4794371.aspx
Leshami
2018/08/07
7690
表分区中的分区交换
插入,更新,删除操作在具有大量数据的表中会变的很慢。通过分区表的分区交换可以快速实现这个过程。
Vincent-yuan
2020/04/08
2.6K0
表分区中的分区交换
浪尖说spark的coalesce的利弊及原理
浪尖的粉丝应该很久没见浪尖发过spark源码解读的文章,今天浪尖在这里给大家分享一篇文章,帮助大家进一步理解rdd如何在spark中被计算的,同时解释一下coalesce降低分区的原理及使用问题。
Spark学习技巧
2020/04/07
4.1K0
Oracle分区表之创建维护分区表索引的详细步骤
墨墨导读:本文来自墨天轮用户投稿,详细描述Oracle分区表之创建维护分区表索引的步骤。
数据和云
2020/07/02
2.4K0
推荐阅读
相关推荐
Greenplum 7 新特性整理
更多 >
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档
本文部分代码块支持一键运行,欢迎体验
本文部分代码块支持一键运行,欢迎体验