samtools軟件從bam文件提取目標(biāo)序列

可以使用samtools fastq命令從bam文件中提取特定的fastq序列葵礼。
使用幫助

samtools fastq
Usage: samtools fastq [options...] <in.bam>
Options:
  -0 FILE              write reads designated READ_OTHER to FILE
  -1 FILE              write reads designated READ1 to FILE
  -2 FILE              write reads designated READ2 to FILE
                       note: if a singleton file is specified with -s, only
                       paired reads will be written to the -1 and -2 files.
  -f INT               only include reads with all  of the FLAGs in INT present [0]
  -F INT               only include reads with none of the FLAGS in INT present [0]
  -G INT               only EXCLUDE reads with all  of the FLAGs in INT present [0]
  -n                   don't append /1 and /2 to the read name
  -N                   always append /1 and /2 to the read name
  -O                   output quality in the OQ tag if present
  -s FILE              write singleton reads designated READ1 or READ2 to FILE
  -t                   copy RG, BC and QT tags to the FASTQ header line
  -T TAGLIST           copy arbitrary tags to the FASTQ header line
  -v INT               default quality score if not given in file [1]
  -i                   add Illumina Casava 1.8 format entry to header (eg 1:N:0:ATCACG)
  -c                   compression level [0..9] to use when creating gz or bgzf fastq files
  --i1 FILE            write first index reads to FILE
  --i2 FILE            write second index reads to FILE
  --barcode-tag TAG    Barcode tag [default: BC]
  --quality-tag TAG    Quality tag [default: QT]
  --index-format STR   How to parse barcode and quality tags

      --input-fmt-option OPT[=VAL]
               Specify a single input file format option in the form
               of OPTION or OPTION=VALUE
      --reference FILE
               Reference sequence FASTA FILE [null]
  -@, --threads INT
               Number of additional threads to use [0]

Reads are designated READ1 if FLAG READ1 is set and READ2 is not set.
Reads are designated READ2 if FLAG READ1 is not set and READ2 is set.
Reads are designated READ_OTHER if FLAGs READ1 and READ2 are either both set
or both unset.
Run 'samtools flags' for more information on flag codes and meanings.

The index-format string describes how to parse the barcode and quality tags, for example:
   i14i8       the first 14 characters are index 1, the next 8 characters are index 2
   n8i14       ignore the first 8 characters, and use the next 14 characters for index 1
If the tag contains a separator, then the numeric part can be replaced with '*' to mean
'read until the separator or end of tag', for example:
   n*i*        ignore the left part of the tag until the separator, then use the second part
               of the tag as index 1

Examples:
 To get just the paired reads in separate files, use:
   samtools fastq -1 paired1.fq -2 paired2.fq -0 /dev/null -s /dev/null -n -F 0x900 in.bam

 To get all non-supplementary/secondary reads in a single file, redirect the output:
   samtools fastq -F 0x900 in.bam > all_reads.fq

Flags對應(yīng)的意義如下:
BAM比對文件中flags對應(yīng)的意義 - 簡書 (jianshu.com)

Flags對應(yīng)的意義.png

以上可以知道數(shù)字4可以將未比對上序列的和比對上序列分開。
數(shù)字64或128可以將read1和read2區(qū)分開的烁。
這是示例序列的比對結(jié)果:

36599658 reads; of these:
36599658 (100.00%) were paired; of these:
36567223 (99.91%) aligned concordantly 0 times
32435 (0.09%) aligned concordantly exactly 1 time
0 (0.00%) aligned concordantly >1 times
----
36567223 pairs aligned concordantly 0 times; of these:
182 (0.00%) aligned discordantly 1 time
----
36567041 pairs aligned 0 times concordantly or discordantly; of these:
73134082 mates make up the pairs; of these:
73133169 (100.00%) aligned 0 times
913 (0.00%) aligned exactly 1 time
0 (0.00%) aligned >1 times
0.09% overall alignment rate

提取所有reads

雙末端測序

ls *.bam |while read id
do
samtools fastq -@ 4 -N $id -1 ${id%%_*}_1.fq -2 ${id%%_*}_2.fq
gzip ${id%%_*}_1.fq ${id%%_*}_2.fq
done

單末端測序

ls *.bam |while read id
do
samtools fastq -@ 2 -N $id > ${id%%_*}.fq
gzip ${id%%_*}.fq
done

如果要提取特定狀態(tài)的reads贴届,需要使用f(提取)或F(過濾)參數(shù)+flags(BAM比對文件中flags對應(yīng)的意義 - 簡書 (jianshu.com)
),例如:-f 2為提取兩條reads都比對上的reads;-F 4:提取所有比對上的reads(包含雙末端中只有單條比對上的reads)盯捌;-F 8過濾兩條reads都沒比對上的reads。
有些狀態(tài)的reads不能通過一個(gè)步驟就將其提饶⒒唷饺著;例如配對的兩條reads只有其中一條比對上滤祖,我想提取另一條∑孔眩可以分為兩步提取:

samtools view -f 4 -b in.bam | samtools fastq -F 8 -1 read1-1.fq -2 read2-2.fq 

或者兩條reads只有其中一條比對上埂材,我想提取比對上的這一條

samtools view -F 4 -b in.bam | samtools fastq -F 2 -1 read1-1.fq -2 read2-2.fq
最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末塑顺,一起剝皮案震驚了整個(gè)濱河市,隨后出現(xiàn)的幾起案子俏险,更是在濱河造成了極大的恐慌严拒,老刑警劉巖,帶你破解...
    沈念sama閱讀 222,681評論 6 517
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件竖独,死亡現(xiàn)場離奇詭異裤唠,居然都是意外死亡,警方通過查閱死者的電腦和手機(jī)莹痢,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 95,205評論 3 399
  • 文/潘曉璐 我一進(jìn)店門种蘸,熙熙樓的掌柜王于貴愁眉苦臉地迎上來,“玉大人竞膳,你說我怎么就攤上這事航瞭。” “怎么了坦辟?”我有些...
    開封第一講書人閱讀 169,421評論 0 362
  • 文/不壞的土叔 我叫張陵刊侯,是天一觀的道長。 經(jīng)常有香客問我锉走,道長滨彻,這世上最難降的妖魔是什么? 我笑而不...
    開封第一講書人閱讀 60,114評論 1 300
  • 正文 為了忘掉前任挪蹭,我火速辦了婚禮亭饵,結(jié)果婚禮上,老公的妹妹穿的比我還像新娘嚣潜。我一直安慰自己冬骚,他們只是感情好,可當(dāng)我...
    茶點(diǎn)故事閱讀 69,116評論 6 398
  • 文/花漫 我一把揭開白布懂算。 她就那樣靜靜地躺著只冻,像睡著了一般。 火紅的嫁衣襯著肌膚如雪计技。 梳的紋絲不亂的頭發(fā)上喜德,一...
    開封第一講書人閱讀 52,713評論 1 312
  • 那天,我揣著相機(jī)與錄音垮媒,去河邊找鬼舍悯。 笑死航棱,一個(gè)胖子當(dāng)著我的面吹牛,可吹牛的內(nèi)容都是我干的萌衬。 我是一名探鬼主播饮醇,決...
    沈念sama閱讀 41,170評論 3 422
  • 文/蒼蘭香墨 我猛地睜開眼,長吁一口氣:“原來是場噩夢啊……” “哼秕豫!你這毒婦竟也來了朴艰?” 一聲冷哼從身側(cè)響起,我...
    開封第一講書人閱讀 40,116評論 0 277
  • 序言:老撾萬榮一對情侶失蹤混移,失蹤者是張志新(化名)和其女友劉穎祠墅,沒想到半個(gè)月后,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體歌径,經(jīng)...
    沈念sama閱讀 46,651評論 1 320
  • 正文 獨(dú)居荒郊野嶺守林人離奇死亡毁嗦,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 38,714評論 3 342
  • 正文 我和宋清朗相戀三年,在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了回铛。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片狗准。...
    茶點(diǎn)故事閱讀 40,865評論 1 353
  • 序言:一個(gè)原本活蹦亂跳的男人離奇死亡,死狀恐怖勺届,靈堂內(nèi)的尸體忽然破棺而出驶俊,到底是詐尸還是另有隱情,我是刑警寧澤免姿,帶...
    沈念sama閱讀 36,527評論 5 351
  • 正文 年R本政府宣布饼酿,位于F島的核電站,受9級特大地震影響胚膊,放射性物質(zhì)發(fā)生泄漏故俐。R本人自食惡果不足惜,卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 42,211評論 3 336
  • 文/蒙蒙 一紊婉、第九天 我趴在偏房一處隱蔽的房頂上張望药版。 院中可真熱鬧,春花似錦喻犁、人聲如沸槽片。這莊子的主人今日做“春日...
    開封第一講書人閱讀 32,699評論 0 25
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽还栓。三九已至,卻和暖如春传轰,著一層夾襖步出監(jiān)牢的瞬間剩盒,已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 33,814評論 1 274
  • 我被黑心中介騙來泰國打工慨蛙, 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留辽聊,地道東北人纪挎。 一個(gè)月前我還...
    沈念sama閱讀 49,299評論 3 379
  • 正文 我出身青樓,卻偏偏與公主長得像跟匆,于是被迫代替她去往敵國和親异袄。 傳聞我的和親對象是個(gè)殘疾皇子,可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 45,870評論 2 361