學(xué)生信的那些事兒之十五 - 生信技能樹sam和bam格式文件的shell小練習(xí)20題

用上一節(jié)下載的bowtie2軟件中自帶的測(cè)試數(shù)據(jù)生成sam和bam文件先:

# 下載bowtie2軟件先
wget -c https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.3.5.1/bowtie2-2.3.5.1-linux-x86_64.zip
unzip bowtie2-2.3.5.1-linux-x86_64.zip
# 找到示例數(shù)據(jù)
cd bowtie2-2.3.5.1-linux-x86_64/example/reads
# 先生成sam文件延届,使用bowtie2命令,絕對(duì)路徑調(diào)用荔燎,參數(shù)什么的DIY吧
/trainee/Jude/practice/bowtie2-2.3.5.1-linux-x86_64/bowtie2 -x /trainee/Jude/practice/bowtie2-2.3.5.1-linux-x86_64/example/index/lambda_virus -1 reads_1.fq -2 reads_2.fq > test.sam
# sam轉(zhuǎn)換成bam
samtools view -bS test.sam > test.bam

一耻姥、統(tǒng)計(jì)共多少條reads(pair-end reads這里算一條)參與了比對(duì)參考基因組

# 這道題的命令很簡(jiǎn)單,但是在寫命令之前花了好長(zhǎng)的時(shí)間學(xué)習(xí)sam文件每一列的含義.下面命令的值是20000有咨,因?yàn)槭请p端琐簇,所以有10000條pair-end reads參與了比對(duì)。
less -SN  test.sam | grep -v  "^@" | wc

二座享、統(tǒng)計(jì)共有多少種比對(duì)的類型(即第二列數(shù)值有多少種)及其分布

# 就是統(tǒng)計(jì)有多少種flag, 取出第二列婉商,排序,然后合并計(jì)數(shù)
less -SN  test.sam | grep -v  "^@" | awk '{print $2}' | sort | uniq -c

三渣叛、篩選出比對(duì)失敗的reads丈秩,看看序列特征

# 先搞明白什么叫作比對(duì)失敗的reads:第6列 ( CIGAR )為*則表示比對(duì)失敗,先計(jì)數(shù)
less -SN  test.sam | grep -v  "^@" | awk '{if($6=="*")print}'| wc

# 再看序列特征, 序列信息在第10例
less -SN  test.sam | grep -v  "^@" | awk '{if($6=="*")print}'| awk '{print$10}' | less -SN

其實(shí)這里有個(gè)疑問,按說第三列如果為*也表示比對(duì)失敶狙谩(沒有比對(duì)上)蘑秽,為什么最終的數(shù)值不一樣呢···

四、比對(duì)失敗的reads區(qū)分成單端失敗和雙端失敗情況箫攀,并且拿到序列ID

# 考的是概念不是技術(shù)啊啊啊
less test.sam |grep -v '^@'|awk '{if($6 =="*") print$1}'|sort -n |uniq -c |grep -w 1
less test.sam |grep -v '^@'|awk '{if($6 =="*") print$1}'|sort -n |uniq -c |grep -w 2

五肠牲、篩選出比對(duì)質(zhì)量值大于30的情況(看第5列)

# 送分題 : 第五列表示比對(duì)質(zhì)量
less test.sam |grep -v '^@'|awk '{if($5 > 30) print}' | less -SN
less test.sam |grep -v '^@'|awk '{if($5 > 30) print}' | wc

六、篩選出比對(duì)成功靴跛,但是并不是完全匹配的序列

# 第一反應(yīng)還是看第6列
less test.sam |grep -v '^@'|awk '{if($6!="*") print$6}'| grep [IDNXSHP]

# 糾結(jié)了很久缀雳,因?yàn)榭傆X得最后的展示結(jié)果應(yīng)該是$10列,畢竟是要看序列嘛汤求。后來一想俏险,單純看序列看不出來這些突變呃。扬绪。竖独。

七、篩選出insert size長(zhǎng)度大于1250bppair-end reads

less test.sam | grep -v '^@'| awk '{if($9>1250) print}'| head | less -SN

八挤牛、統(tǒng)計(jì)參考基因組上面各條染色體的成功比對(duì)reads數(shù)量

# 成功比對(duì)后第3列會(huì)顯示對(duì)應(yīng)的參考序列信息莹痢,無法比對(duì)則顯示*,所以取出第3列,排序然后計(jì)數(shù)
less test.sam | grep -v '^@'| cut -f 3 | sort -u | uniq -c

九竞膳、篩選出原始fq序列里面有N的比對(duì)情況

# 比對(duì)情況看第6列航瞭,不過首先需要篩選出有N的行
less test.sam | grep -v '^@'| awk '{if($10~N) print$6}' | less

十、篩選出原始fq序列里面有N坦辟,但是比對(duì)的時(shí)候卻是完全匹配的情況

# 同上題思路
less test.sam | grep -v '^@'| awk '{if($10~N) print}' | awk '{if($6!~"[IDNXSHP*]") print}' | less -SN

發(fā)現(xiàn)卡殼主要是因?yàn)椴皇煜wk的語法和SAM數(shù)據(jù)結(jié)構(gòu)刊侯,這樣的時(shí)候千萬不要跟自己較勁,去看別人的答案和教程吧锉走,因?yàn)椴皇沁壿嬐评淼膯栴}滨彻,而更像是公式,參考別人的挪蹭,學(xué)會(huì)了就好亭饵。

十一、sam文件里面的頭文件行數(shù)

# 頭文件結(jié)構(gòu)是以@符號(hào)打頭
less test.sam | grep '^@' | wc

十二梁厉、sam文件里每一行的tags個(gè)數(shù)一樣嗎

# 搞清楚tags是啥辜羊,在什么地方先:在第12列及后面的列,為可選字段(optional fields)词顾,格式如:TAG:TYPE:VALUE八秃,其中TAG有兩個(gè)大寫字母組成,每個(gè)TAG代表一類信息计技,每一行一個(gè)TAG只能出現(xiàn)一次喜德,TYPE表示TAG對(duì)應(yīng)值的類型,可以是字符串垮媒、整數(shù)、字節(jié)航棱、數(shù)組等睡雇。

less test.sam | grep -v '^@' | cut -f 12- | less -SN

# 肉眼可見的不一樣

十三、sam文件里每一行的tags個(gè)數(shù)分別是多少個(gè)

# 跟上題一樣, 由輸出結(jié)果可知絕大部分為9個(gè)饮醇,也有2個(gè)它抱,8個(gè)的
less test.sam | grep -v '^@' | cut -f 12- | awk -F '\t' '{print "columns:" NF }'

這道題做出來很開心,因?yàn)槟壳澳芩训降钠渌藢懙拇鸢高@道題都是空著的朴艰,哈哈哈~~

十四观蓄、sam文件里記錄的參考基因組染色體長(zhǎng)度分別是?

# 參考基因組染色體的長(zhǎng)度信息記錄在哪里祠墅?
less test.sam | grep '^@' | grep LN

十五侮穿、找到比對(duì)情況有insertion情況的

# 有insertion表示第6列里面有I
less test.sam | grep -v '^@'| cut -f 6 | grep I 

十六、找到比對(duì)情況有deletion情況的

# 有deletion表示第6列里面有D
less test.sam | grep -v '^@'| cut -f 6 | grep D

十七毁嗦、取出位于參考基因組某區(qū)域的比對(duì)記錄亲茅,比如 5013到50130 區(qū)域

# 比對(duì)區(qū)域看第4列,表示比對(duì)的起始位置
less test.sam | grep -v '^@'| awk '{if($4>5013 && $4<50130) print}' | less -SN

十八、把sam文件按照染色體以及起始坐標(biāo)排序

# 就看sort怎么用克锣,首先以染色體排序:
 less test.sam | grep -v '^@' | sort -k 3 | less -S
 
# 以其實(shí)坐標(biāo)排序茵肃,坐標(biāo)是數(shù)字,需-n參數(shù)
 less test.sam | grep -v '^@' | sort -n -k 4 | less -S

查看結(jié)果可知如果按照升序排列的話最前面的都是星號(hào)(*)袭祟,可以改為降序验残。

十九、找到 102M3D11M 的比對(duì)情況巾乳,計(jì)算其reads片段長(zhǎng)度

# grep一下您没,然后cut或者awk一下
less test.sam | grep "102M3D11M" | cut -f 10 | wc

二十、安裝samtools軟件后使用samtools軟件的各個(gè)功能嘗試把上述題目重新做一遍

其它概念題

  1. 高通量測(cè)序數(shù)據(jù)分析中想鹰,序列與參考序列的比對(duì)后得到的標(biāo)準(zhǔn)文件為什么文件紊婉?

? SAM文件

  1. sam文件是一種比對(duì)后的文件,能直接查看嗎?

? 可以辑舷,用cat,less等命令均可查看

  1. Sam/Bam文件分為兩部分喻犁,一部分以@開頭的是什么序列,另一部分是什么序列何缓?

? 以@開頭的是頭部區(qū)域肢础,記錄比對(duì)的整體信息;另一部分為主體區(qū)碌廓,記錄詳細(xì)的比對(duì)信息传轰,至少11列

  1. Sam文件除頭文件,以什么符號(hào)分割文本的谷婆,比對(duì)部分信息一行是多少列慨蛙?你能用awk計(jì)算列數(shù)嗎?

? 分隔符為tab纪挎,比對(duì)部分的列數(shù)必須有的是11列期贫,剩余的列數(shù)不一定;用awk計(jì)算列數(shù)參考第13題答案

  1. Sam/Bam文件的@SQ開頭的行是什么异袄?你能生成一個(gè)文本通砍,兩列,一列是參考基因組染色體/sca id烤蜕,一列是長(zhǎng)度嗎封孙?

? 是參考基因組序列信息

  1. Sam文件的比對(duì)位置是從1還是0開始的?

? 從1開始

  1. 常見CIGAR 字符串各字母代表的意義

? “M”表示 match或 mismatch讽营;

? “I”表示 insert

? “D”表示 deletion

? “N”表示 skipped(跳過這段區(qū)域)

? “S”表示 soft clipping(被剪切的序列存在于序列中)

? “H”表示 hard clipping(被剪切的序列不存在于序列中)

? “P”表示 padding

? “=”表示 match

? “X”表示 mismatch(錯(cuò)配虎忌,位置是一一對(duì)應(yīng)的)

  1. 比對(duì)質(zhì)量的數(shù)字是哪一列?越大比對(duì)質(zhì)量越好還是越小越好斑匪?

? 第五列呐籽,越大越好锋勺,60最好

  1. Sam文件的flag是第幾列,數(shù)字代表什么意義狡蝶?是怎么計(jì)算來的庶橱?

? 第二列,代表比對(duì)的情況贪惹,具體計(jì)算看表吧

  1. Sam文件怎么轉(zhuǎn)bam文件苏章?用什么指令?為什么要轉(zhuǎn)換?

? 命令參考samtools view -bS test.sam > test.bam奏瞬;Bam文件相對(duì)于Sam文件進(jìn)一步壓縮枫绅,同樣的空間可以儲(chǔ)存更多數(shù)據(jù)

  1. Bam文件查看用什么指令?如果需要查看頭文件需要增加參數(shù)硼端?

? 指令: samtools view; 查看頭文件指令:samtools view test.bam | grep "^@"

  1. Bam文件為什么要排序并淋?排序后的bam和未排序的bam頭文件和chr position列有什么區(qū)別?

? 這個(gè)不是很確定···

  1. Bam文件建索引的指令是什么珍昨?

? samtools index test.bam test.bam.bai

  1. Bam文件可視化用什么工具县耽?查看時(shí)需要建立索引嗎?

? 可視化可以用IGV吧

  1. Bam文件統(tǒng)計(jì)reads比對(duì)情況用samtools的哪一個(gè)子命令镣典?

? samtools tview

  1. SE測(cè)序和PE測(cè)序的所比對(duì)后得到的sam文件的區(qū)別在哪里兔毙?

? 根據(jù)對(duì)SAM文件格式的理解,F(xiàn)lag值會(huì)不一樣兄春;其他的不太確定···

  1. Bam能用gzip再壓縮嗎澎剥?

? not sure about this

  1. Sam文件通常由哪些比對(duì)軟件得到的?

? Bowtie2啊赶舆,Hisat啊什么的

  1. Sam/Bam文件能轉(zhuǎn)成fastq文件嗎哑姚?

? 當(dāng)然可以

  1. 有時(shí)候不能通過文件名的后綴來區(qū)別是sam還是bam逆瑞,你是怎么區(qū)別的纯趋?

? Sam文件可以直接用cat蓬抄,less命令查看陷遮,Bam為二進(jìn)制文件,需要用samtools view命令查看

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
  • 序言:七十年代末抄课,一起剝皮案震驚了整個(gè)濱河市,隨后出現(xiàn)的幾起案子,更是在濱河造成了極大的恐慌蒸辆,老刑警劉巖,帶你破解...
    沈念sama閱讀 217,277評(píng)論 6 503
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件析既,死亡現(xiàn)場(chǎng)離奇詭異躬贡,居然都是意外死亡,警方通過查閱死者的電腦和手機(jī)眼坏,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 92,689評(píng)論 3 393
  • 文/潘曉璐 我一進(jìn)店門拂玻,熙熙樓的掌柜王于貴愁眉苦臉地迎上來,“玉大人,你說我怎么就攤上這事檐蚜∑嵌” “怎么了?”我有些...
    開封第一講書人閱讀 163,624評(píng)論 0 353
  • 文/不壞的土叔 我叫張陵闯第,是天一觀的道長(zhǎng)市栗。 經(jīng)常有香客問我,道長(zhǎng)咳短,這世上最難降的妖魔是什么填帽? 我笑而不...
    開封第一講書人閱讀 58,356評(píng)論 1 293
  • 正文 為了忘掉前任,我火速辦了婚禮咙好,結(jié)果婚禮上篡腌,老公的妹妹穿的比我還像新娘。我一直安慰自己勾效,他們只是感情好嘹悼,可當(dāng)我...
    茶點(diǎn)故事閱讀 67,402評(píng)論 6 392
  • 文/花漫 我一把揭開白布。 她就那樣靜靜地躺著葵第,像睡著了一般绘迁。 火紅的嫁衣襯著肌膚如雪。 梳的紋絲不亂的頭發(fā)上卒密,一...
    開封第一講書人閱讀 51,292評(píng)論 1 301
  • 那天缀台,我揣著相機(jī)與錄音,去河邊找鬼哮奇。 笑死膛腐,一個(gè)胖子當(dāng)著我的面吹牛,可吹牛的內(nèi)容都是我干的鼎俘。 我是一名探鬼主播哲身,決...
    沈念sama閱讀 40,135評(píng)論 3 418
  • 文/蒼蘭香墨 我猛地睜開眼,長(zhǎng)吁一口氣:“原來是場(chǎng)噩夢(mèng)啊……” “哼贸伐!你這毒婦竟也來了勘天?” 一聲冷哼從身側(cè)響起,我...
    開封第一講書人閱讀 38,992評(píng)論 0 275
  • 序言:老撾萬榮一對(duì)情侶失蹤捉邢,失蹤者是張志新(化名)和其女友劉穎脯丝,沒想到半個(gè)月后,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體伏伐,經(jīng)...
    沈念sama閱讀 45,429評(píng)論 1 314
  • 正文 獨(dú)居荒郊野嶺守林人離奇死亡宠进,尸身上長(zhǎng)有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 37,636評(píng)論 3 334
  • 正文 我和宋清朗相戀三年,在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了藐翎。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片材蹬。...
    茶點(diǎn)故事閱讀 39,785評(píng)論 1 348
  • 序言:一個(gè)原本活蹦亂跳的男人離奇死亡实幕,死狀恐怖,靈堂內(nèi)的尸體忽然破棺而出堤器,到底是詐尸還是另有隱情昆庇,我是刑警寧澤,帶...
    沈念sama閱讀 35,492評(píng)論 5 345
  • 正文 年R本政府宣布吼旧,位于F島的核電站凰锡,受9級(jí)特大地震影響,放射性物質(zhì)發(fā)生泄漏圈暗。R本人自食惡果不足惜掂为,卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 41,092評(píng)論 3 328
  • 文/蒙蒙 一、第九天 我趴在偏房一處隱蔽的房頂上張望员串。 院中可真熱鬧勇哗,春花似錦、人聲如沸寸齐。這莊子的主人今日做“春日...
    開封第一講書人閱讀 31,723評(píng)論 0 22
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽渺鹦。三九已至扰法,卻和暖如春,著一層夾襖步出監(jiān)牢的瞬間毅厚,已是汗流浹背塞颁。 一陣腳步聲響...
    開封第一講書人閱讀 32,858評(píng)論 1 269
  • 我被黑心中介騙來泰國(guó)打工, 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留吸耿,地道東北人祠锣。 一個(gè)月前我還...
    沈念sama閱讀 47,891評(píng)論 2 370
  • 正文 我出身青樓,卻偏偏與公主長(zhǎng)得像咽安,于是被迫代替她去往敵國(guó)和親伴网。 傳聞我的和親對(duì)象是個(gè)殘疾皇子,可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 44,713評(píng)論 2 354

推薦閱讀更多精彩內(nèi)容