生物信息分析:从入门到精通生信笔记比对

bwa samtools bowtie2大杂烩

2018-09-27  本文已影响53人  11的雾

比对,先要建index

bwa index /cygene/database/TCR/TRAVBV.seq.fa
一般bwa比对完后需要用samtools转换成bam或者排序什么的,都需要用samtools首先对ref进行建index
samtools faidx TRAVBV.seq.fa
基础参数介绍

-R STR        read group header line such as '@RG\tID:foo\tSM:bar'  [null]

如: "@RG\tID:0\tPL:BGISEQ-500\tPU:0\tLB:UN\tSM:0\tCN:BGI"

bwa mem 搭配samtools比对,及sort排序:

PE(paired end)

bwa mem -t 18 /cygene/database/TCR/TCR.all.vdjab.fa /cygene/data/ANCCR180494_PM-CR180494-04_AHNTL2CCXY_2018-10-01/Cleandata/RP01G9E1L3/RP01G9E1L3_R1.fq.gz /cygene/data/ANCCR180494_PM-CR180494-04_AHNTL2CCXY_2018-10-01/Cleandata/RP01G9E1L3/RP01G9E1L3_R2.fq.gz |samtools view -buhS -t /cygene/database/TCR/TCR.all.vdjab.fa.fai - |samtools sort -n -o Cleandata_sample1_Rd12.vdjab.bam -

SE (single end)

bwa mem -t 18 /cygene/database/TCR/TCR.all.vdjab.fa /cygene/data/ANCCR180494_PM-CR180494-03_000000000-C2R58_2018-09-14/Rawdata/RP01G9E1L1/RP01G9E1L1_R1.fq.gz |samtools view -buhS -t /cygene/database/TCR/TCR.all.vdjab.fa.fai - |samtools sort -n -o sample1_R1.vdjab.bam -

bwa aln

bwa aln用于短read,
PE:比对两次,再合并
SE:

比对结果处理

比对结果中没有比对上的flag为4
查看bam/sam文件用 samtools view xxx.bam |le
过滤没有比对上的结果 用-F:

-F INT   only include reads with none of the FLAGS in INT present [0]

例子1:从bam中过滤掉没有比对上的信息,并将比对上的部分保存到sam中,
samtools view -F 4 xxx.bam >xxx.mapped.sam
例子2:从bam中过滤没有比对上的信息,并保存到bam中:(要加-h,表示输出header)
samtools view -F 4 -h xxx.bam |samtools view -h -o xxx.mapped.bam -
注意:虽然有些没有比对上的结果 的flag值是141,77等,但是都可以用4过滤掉。

image.png

========================================
提取为fastq:主要参数为-f,-F,-G三个。
提取没有比对上的信息到fastq用:
samtools fastq -f 4 xxx.bam > xxx.unmapped.fastq
提取比对上的reads到fastq用
samtools fastq -F 4 xxx.bam > xxx.mapped.fastq
查询flag值的意义:https://broadinstitute.github.io/picard/explain-flags.html

samtools sort 报错案例:

使用samtools 为bam建立index之前,要对bam进行sort,而且sort时候不能加-n参数(sort by read name)。否则samtools index会报错:

[E::hts_idx_push] NO_COOR reads not in a single block at the end 280 -1
samtools index: failed to create index for "sample1_R1.vdjab.bam"

bam与sam格式相互转换

BAM 转SAM :samtools view -h -o out.sam out.bam
SAM转BAM :samtools view -bS out.sam >out.bam
-b 意思是输出为BAM format
-S 意思是输入为SAM,如果@SQ 缺剩, 要写-t,所以如果没有@SQ
samtools faidx ref.fa
samtools view -bt ref.fa.fai out.sam > out.bam

bowtie2

bowtie2建index用:
bowtie2-build -f TRAVBV.seq.fa TRAVBV.seq前缀就可以
全局局部比对要加参数--local
只比对正链或只比对负链用--norc,--nofw
single end数据用-U比对实例:
bowtie2 -p 20 -x /cygene/database/TCR/TRAVBV.seq -U /cygene/data/ANCCR180494_PM-CR180494-03_000000000-C2R58_2018-09-14/Rawdata/RP01G9E1L1/RP01G9E1L1_R1.fq.gz -S sample1_r1_bowtie.fw.sam --nofw --local

上一篇下一篇

猜你喜欢

热点阅读