生信小工具专题:BBTools/BBMap Suite 的使用 (3)

接着上一次内容继续给大家带来实用的BBTools的小教程。

BBMap suite - Simulation of data (数据模拟)

除了之前介绍的工具,BBMap还可以用来,
从参考基因组产生随机合reads。reads的名称表示其基因组来源。并可以允许准确定的insert size和合成突变的类型,大小和速率等具体内容。 Illumina和PacBio类型的reads都可以使用randomreads这个小工具来生成。

您可以指定PE 的reads,insert size的分布,reads长度(或长度范围)等。 Randomreads专门设计用于对突变进行出色的控制。你可以准确地或概率地指定每个reads所包含的snps,insertion,deletions和Ns的数量;这些突变的长度可以进行单独的定制,质量值可以交替设置,以允许根据质量生成错误。所有reads都使用其基因组来源进行注释。

请记住,50%的reads将来自正链,50%来自负链。因此,reads将完全匹配参考基因组,或者它的反向补链将完美匹配。

基本用法

randomreads.sh ref=<file> out=<file> length=<number> reads=<number>

介绍那么多,可能大家还是会有点疑惑究竟,这个工具有啥用?废话不多说,通过几行代码带你快速了解该工具的用处:

  1. 模拟高通量测序Illumina数据集。

首先先下载序列号为NC_010468 fasta格式的文件

efetch -db=nuccore -format=fasta -id=NC_010468 > NC_010468.fa

我们将通过以下命令简化NC_010468.fa中的fasta标头,以便新标头读取> NC_010468_E_coli

sed '1 s/^.*$/\>NC_010468_E_coli/' NC_010468.fa > NC_010468_new.fa

我们将使用此参考序列文件生成的合成illumina配对末端数据集(2 x 300 bp)

randomreads.sh ref=NC_010468_new.fa out=NC_010468_synth.fq.gz length=300 paired=t mininsert=100 maxinsert=400 reads=1000000
 
#默认情况下,生成的reads采用交错格式。它们可以通过`reformat.sh`进一步分成单独的read文件

reformat.sh in=NC_010468_synth.fq.gz out1=NC_010468_synth_R1.fq.gz out2=NC_010468_synth_R2.fq.gz addcolon
  1. 使用randomreads生成PacBio格式数据集(使用PacBio错误模型)。

模拟一个包含10,000个读数的PacBio数据集,其最小读取长度为1kb,最大读取长度为5kb。

randomreads.sh ref=NC_010468_new.fa out=NC_010468_synth_pacbio.fq.gz minlength=1000 maxlength=5000 reads=10000 pacbio=t
  1. 生成更复杂/有趣的合成数据集,可以模拟宏基因组。
    下载数据
efetch -db=nuccore -format=fasta -id=NC_009614 > NC_009614.fa
efetch -db=nuccore -format=fasta -id=NC_013520 > NC_013520.fa

metagenome将为每个scaffold分配一个随机指数变量,该变量决定从该scaffold生成read的概率。因此,如果将20个细菌基因组连接在一起,则可以运行随机read并获得类似宏基因组的分布。当使用转录组参考时,它也可以用于RNA-seq。

覆盖率取决于每个参考序列水平,因此如果细菌组装具有多个contigs,您可能需要先将它们与fuse.sh粘在一起,然后再与其他参考物连接。

首先将ref 合成单一的文件:

cat NC_010468.fa NC_013520.fa NC_009614.fa AF086833.fa > metag.fa

然后我们将使用它作为宏基因组生成的输入,然后使用reformat.sh将读取分成两个文件。

randomreads.sh ref=metag.fa out=stdout length=300 paired=t mininsert=100 maxinsert=400 reads=5000000 metagenome=t coverage=3 | reformat.sh in=stdin out1=metag_R1.fq.gz out2=metag_R2.fq.gz interleaved addslash

重新格式化和过滤序列数据

Reformat.sh这个工具包前面一直有用到,接下来会详细讲解一下他的使用。

功能:
格式转换/操作/分析序列数据。

重新格式化reads以更改ASCII质量编码,交错,文件格式或压缩格式。可选择执行其他功能,如质量修整,子集和子采样。支持fastq,fasta,fasta,bam,gzip 等得格式。

基本用法:

reformat.sh in=<file> in2=<file2> out=<outfile> out2=<outfile2>
#in2 and out2 这两个options是PE reads的文件才会用到的(双末端)

具体应用

  1. 切除一部分的reads文件

使用reformat.sh可以轻松地从文件中产生随机获取的reads。

#随机抽取300个reads
reformat.sh in=SRR1972739_1.fastq out=SRR1972739_1_sample.fastq sample=300


#您可以通过执行以总reads的百分比进行采样(以下示例中为10%):
reformat.sh in=SRR1972739_1.fastq out=SRR1972739_1_10perc.fastq samplerate=0.1
  1. 验证配对端数据的reads是正确配对的。
reformat.sh in1=r1.fq.gz in2=r2.fq.gz vpair
  1. 按reads的长度过滤SAM / BAM文件
reformat.sh in=x.sam out=y.sam minlength=50 maxlength=200
  1. 过滤/检测q SAM / BAM文件中的拼接reads

您可以将maxdellen(在下面的示例中为50)设置为你认为最小的任何长度的eletion事件,以表示拼接,这取决于你的物种的类型。

reformat.sh in=mapped.bam out=filtered.bam maxdellen=50

从SAM / BAM文件中获取比对不上的reads

reformat.sh in=all.sam out=unmapped.fq.gz unmappedonly
repair.sh in=unmapped.fq.gz  out=r1.fq.gz out2=r2.fq.gz outs=singleton.fq

删除重复的序列数据集

从一组文件中删除重复的读取/序列。 类似picard tool的Markduplicates的功能。

dedupe.sh in=<file or stdin> out=<file or stdout>

在运行完成后,还会生成一组全面的统计信息。它们通常被写入\但可以很容易地捕获到日志文件中(例如下面的例子)

Input:                          15000 reads             1515000 bases.
Duplicates:                     7445 reads (49.63%)     751945 bases (49.63%)           0 collisions.
Containments:                   0 reads (0.00%)         0 bases (0.00%)         113772 collisions.
Result:                         7555 reads (50.37%)     763055 bases (50.37%)
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,222评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,455评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,720评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,568评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,696评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,879评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,028评论 3 409
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,773评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,220评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,550评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,697评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,360评论 4 332
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,002评论 3 315
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,782评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,010评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,433评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,587评论 2 350

推荐阅读更多精彩内容