10.12 bwa使用 安装文件路径与使用 sh权限

我们这里将用于流程构建的BWA就是其中最优秀的一个,它将BW(Burrows-Wheeler)压缩算法和后缀树相结合,能够让我们以较小的时间和空间代价,获得准确的序列比对结果。

别人的已安装文件打包传递后使用:

1、连接服务器

2、家目录下,Users,ls -all (或者 打开/etc里profile)

3、vim .bashrc

```

export PATH="your path:$PATH"

```

添加的路径是绿色的可执行文件所在的位置

4、source .bashrc 

5、在全局中输入bwa测试是否可以使用

自己安装的软件路径都已经配置好了,所以不需要这样操作就可以直接使用。


BWA MEM比对模块是有一定适用范围的:它是专门为长read比对设计的,目的是为了解决,第三代测序技术这种能够产生长达几十kb甚至几Mbp的read情况。一般只有当read长度≥70bp的时候,才推荐使用,如果比这个要小,建议使用BWA ALN模块。

BWA ALN模块

1: 建立 Index

根据reference genome data(e.g. reference.fa) 建立 Index File

    bwa index -a bwtsw reference.fa

(bwa index 可以查看更多使用方法)

step 2: 寻找 SA coordinates

如果是pair-end 数据(leftRead.fastq和rightRead.fastq)两个文件分别处理

    bwa aln reference.fa leftRead.fastq > leftRead.sai

    bwa aln reference.fa rightRead.fastq > rightRead.sai

    bwa aln reference.fa singleRead.fastq > singleRead.sai

如果希望多线程运行,在其中加入 -t这个参数,另外-f这个参数可以指定结果输出文件,如:

    bwa aln -c -t 3 -f leftreads.sai reference.fa leftreads.fastq

step 3:转换SA coordinates输出为sam

如果是pair-end数据

    bwa sampe -f pair-end.sam reference.fa leftRead.sai rightRead.sai leftRead.fastq rightread.fastq

如果是single reads数据

    bwa samse -f single.sam reference.fa single.sai single.fastq

bwa mem模块

https://www.jianshu.com/p/859c0345624c

bwa mem [options] <idxbase> <in1.fq> [in2.fq]

bwa mem -t 4 -R'@RG\tID:foo_lane\tPL:illumina\tLB:library\tSM:sample_name'/path/to/human.fasta read_1.fq.gz read_2.fq.gz > sample_name.sam

-t,线程数,我们在这里使用4个线程;-R 接的是Read Group的字符串信息,这是一个非常重要的信息,以@RG开头,它是用来将比对的read进行分组的。不同的组之间测序过程被认为是相互独立的,这个信息对于我们后续对比对数据进行错误率分析和Mark duplicate时非常重要。在Read Group中,有如下几个信息非常重要:

(1) ID,这是Read Group的分组ID,一般设置为测序的lane ID(不同lane之间的测序过程认为是独立的),下机数据中我们都能看到这个信息的,一般都是包含在fastq的文件名中;

(2) PL,指的是所用的测序平台,这个信息不要随便写!特别是当我们需要使用GATK进行后续分析的时候,更是如此!这是一个很多新手都容易忽视的一个地方,在GATK中,PL只允许被设置为:ILLUMINA,SLX,SOLEXA,SOLID,454,LS454,COMPLETE,PACBIO,IONTORRENT,CAPILLARY,HELICOS或UNKNOWN这几个信息。基本上就是目前市场上存在着的测序平台,当然,如果实在不知道,那么必须设置为UNKNOWN,名字方面不区分大小写。如果你在分析的时候这里没设置正确,那么在后续使用GATK过程中可能会碰到类似如下的错误。

(3) SM,样本ID,同样非常重要,有时候我们测序的数据比较多的时候,那么可能会分成多个不同的lane分布测出来,这个时候SM名字就是可以用于区分这些样本。

(4) LB,测序文库的名字,这个重要性稍微低一些,主要也是为了协助区分不同的group而存在。文库名字一般可以在下机的fq文件名中找到,如果上面的lane ID足够用于区分的话,也可以不用设置LB;

除了以上这四个之外,还可以自定义添加其他的信息,不过如无特殊的需要,对于序列比对而言,这4个就足够了。这些信息设置好之后,在RG字符串中要用制表符(\t)将它们分开

sam转换bam(sam特殊二进制)

bwa mem -t 4 -R'@RG\tID:foo_lane\tPL:illumina\tLB:library\tSM:sample_name'/path/to/human.fasta read_1.fq.gz read_2.fq.gz | samtools view -S -b - > sample_name.bam



time$bwamem -t 40 -M -Y -R"@RG\tID:$RGID\tPL:ILLUMINA\tSM:$sample"$reference\$outdir/cleanfq/${fq_file_name}.paired.1.fq.gz$outdir/cleanfq/${fq_file_name}.paired.2.fq.gz |$samtoolsview -Sb - >$outdir/bwa/${sample}.bam && \echo"** BWA MEM done **"&& \

-M 将 shorter split hits 标记为次优,以兼容 Picard’s markDuplicates 软件。

-Y Use soft clipping CIGAR operation for supplementary alignments. By default, BWA-MEM uses soft clipping for the primary alignment and hard clipping for supplementary alignments.

链接:https://www.jianshu.com/p/d6cc153be28d



设置权限 chmod 777 ./*.sh

运行


read比对结束,成功得bam文件

以上,我们就完成了read比对的步骤。接下来是排序:


排序


-@,用于设定排序时的线程数;-m,每个线程排序时最大的内存消耗;-O 指定输出为bam格式;-o 是输出文件的名字.

在排序完成之后我们就可以开始执行去除重复(准确来说是 去除PCR重复序列)的步骤了。


done

把参数REMOVE_DUPLICATES设置为ture,那么重复序列就被删除掉,不会在结果文件中留存。我比较建议使用第一种做法,只是标记出来,并留存这些序列,以便在你需要的时候还可以对其做分析。

这一步完成之后,我们需要为sample_name.sorted.markdup.bam创建索引文件,它的作用能够让我们可以随机访问这个文件中的任意位置,而且后面的“局部重比对”步骤也要求这个BAM文件一定要有索引,命令如下:

$ samtools index sample_name.sorted.markdup.bam

merge

接下来是局部区域重比对,通常也叫Indel局部区域重比对。有时在进行这一步骤之前还有一个merge的操作,将同个样本的所有比对结果合并成唯一一个大的BAM文件【注】,merge的例子如下:

$ samtools merge <out.bam> <in1.bam> [<in2.bam> ... <inN.bam>]

【注意】之所以会有这种情况,是因为有些样本测得非常深,其测序结果需要经过多次测序(或者分布在多个不同的测序lane中)才全部获得,这个时候我们一般会先分别进行比对并去除重复序列后再使用samtools进行合并。

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

推荐阅读更多精彩内容

  • wes定义: 全外显子组测序,是利用目标序列捕获技术, 将全基因组编码基因外显子区域的DNA捕获并富集后,进行高通...
    凤凰_0949阅读 4,144评论 0 7
  • 在上一篇文章《使用MutMap快速定位突变基因:原理篇》中,我对MutMap的原理进行的简单的介绍。在本篇教程中,...
    曹务强阅读 4,438评论 1 15
  • SAM分为两部分,注释信息(header section)和比对结果部分(alignment section)。 ...
    苏牧传媒阅读 7,504评论 0 2
  • 有一瞬间不想活了。 我告诉我朋友说“假如有一天,我找不到生活的乐趣了,只是为了钱而存活着的时候我就去自杀。” 朋友...
    林适情阅读 840评论 0 0
  • 那天早上起床,妈妈问我早上吃什么,我说:“我想去外面买早点吃。”于是妈妈答应了我的要求。 爸爸骑着...
    阳光成子阅读 311评论 0 6