用上一節(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)度大于1250bp的 pair-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è)功能嘗試把上述題目重新做一遍
其它概念題
- 高通量測(cè)序數(shù)據(jù)分析中想鹰,序列與參考序列的比對(duì)后得到的標(biāo)準(zhǔn)文件為什么文件紊婉?
? SAM文件
- sam文件是一種比對(duì)后的文件,能直接查看嗎?
? 可以辑舷,用cat,less等命令均可查看
- Sam/Bam文件分為兩部分喻犁,一部分以@開頭的是什么序列,另一部分是什么序列何缓?
? 以@開頭的是頭部區(qū)域肢础,記錄比對(duì)的整體信息;另一部分為主體區(qū)碌廓,記錄詳細(xì)的比對(duì)信息传轰,至少11列
- Sam文件除頭文件,以什么符號(hào)分割文本的谷婆,比對(duì)部分信息一行是多少列慨蛙?你能用awk計(jì)算列數(shù)嗎?
? 分隔符為tab纪挎,比對(duì)部分的列數(shù)必須有的是11列期贫,剩余的列數(shù)不一定;用awk計(jì)算列數(shù)參考第13題答案
- Sam/Bam文件的@SQ開頭的行是什么异袄?你能生成一個(gè)文本通砍,兩列,一列是參考基因組染色體/sca id烤蜕,一列是長(zhǎng)度嗎封孙?
? 是參考基因組序列信息
- Sam文件的比對(duì)位置是從1還是0開始的?
? 從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)的)
- 比對(duì)質(zhì)量的數(shù)字是哪一列?越大比對(duì)質(zhì)量越好還是越小越好斑匪?
? 第五列呐籽,越大越好锋勺,60最好
- Sam文件的flag是第幾列,數(shù)字代表什么意義狡蝶?是怎么計(jì)算來的庶橱?
? 第二列,代表比對(duì)的情況贪惹,具體計(jì)算看表吧
- Sam文件怎么轉(zhuǎn)bam文件苏章?用什么指令?為什么要轉(zhuǎn)換?
? 命令參考samtools view -bS test.sam > test.bam奏瞬;Bam文件相對(duì)于Sam文件進(jìn)一步壓縮枫绅,同樣的空間可以儲(chǔ)存更多數(shù)據(jù)
- Bam文件查看用什么指令?如果需要查看頭文件需要增加參數(shù)硼端?
? 指令: samtools view; 查看頭文件指令:samtools view test.bam | grep "^@"
- Bam文件為什么要排序并淋?排序后的bam和未排序的bam頭文件和chr position列有什么區(qū)別?
? 這個(gè)不是很確定···
- Bam文件建索引的指令是什么珍昨?
? samtools index test.bam test.bam.bai
- Bam文件可視化用什么工具县耽?查看時(shí)需要建立索引嗎?
? 可視化可以用IGV吧
- Bam文件統(tǒng)計(jì)reads比對(duì)情況用samtools的哪一個(gè)子命令镣典?
? samtools tview
- SE測(cè)序和PE測(cè)序的所比對(duì)后得到的sam文件的區(qū)別在哪里兔毙?
? 根據(jù)對(duì)SAM文件格式的理解,F(xiàn)lag值會(huì)不一樣兄春;其他的不太確定···
- Bam能用gzip再壓縮嗎澎剥?
? not sure about this
- Sam文件通常由哪些比對(duì)軟件得到的?
? Bowtie2啊赶舆,Hisat啊什么的
- Sam/Bam文件能轉(zhuǎn)成fastq文件嗎哑姚?
? 當(dāng)然可以
- 有時(shí)候不能通過文件名的后綴來區(qū)別是sam還是bam逆瑞,你是怎么區(qū)別的纯趋?
? Sam文件可以直接用cat蓬抄,less命令查看陷遮,Bam為二進(jìn)制文件,需要用samtools view命令查看