【生信技能樹】sam和bam格式文件的shell小練習(xí)
LINUX練習(xí)題
1) 統(tǒng)計共多少條reads(pair-end reads這里算一條)參與了比對參考基因組
$grep -v '^@' tmp.sam |cut -f 1|sort -V | uniq -c | wc -l
10000
- 統(tǒng)計共有多少種比對的類型(即第二列數(shù)值有多少種)及其分布。
$grep -v '^@' tmp.sam | cut -f 2 | sort -V |uniq -c
24 65
165 69
153 73
213 77
2 81
4650 83
136 89
4516 99
125 101
16 113
24 129
153 133
165 137
213 141
4516 147
125 153
2 161
4650 163
136 165
16 177
$grep -v '^@' tmp.sam |cut -f 2|sort -V | uniq -c | wc -l
20
3）篩選出比對失敗的reads芙盘，看看序列特征谤逼。
$grep -v '^@' tmp.sam | awk 'match($3,"*"){print;}'
$grep -v '^@' tmp.sam | awk 'match($3,"*"){print;}' | wc -l
426
$grep -v '^@' tmp.sam | awk 'match($3,"*"){print;}' | head -10 |less -SN
- 比對失敗的reads區(qū)分成單端失敗和雙端失敗情況豫缨，并且拿到序列ID
$grep -v '^@' tmp.sam | awk 'match($3,"*"){print($1);}' |\
sort -V | uniq -c | awk '$1==2{print($2);}' #雙端失敗
$grep -v '^@' tmp.sam | awk 'match($3,"*"){print($1);}' |\
sort -V | uniq -c | awk '$1==2{print($2);}' | wc -l
213
$grep -v '^@' tmp.sam | awk 'match($3,"*"){print($1);}' |\
sort -V | uniq -c | awk '$1==1{print($2);}' #單端失敗
$grep -v '^@' tmp.sam | awk 'match($3,"*"){print($1);}' |\
sort -V | uniq -c | awk '$1==1{print($2);}' | wc -l
0
- 篩選出比對質(zhì)量值大于30的情況（看第5列）
$grep -v '^@' tmp.sam | awk '$5>30{print;}'
$grep -v '^@' tmp.sam | awk '$5>30{print;}' | wc -l
18632
$grep -v '^@' tmp.sam | awk '$5>30{print;}'| head -10 |less -SN
- 篩選出比對成功，但是并不是完全匹配的序列
- 篩選出inset size長度大于1250bp的 pair-end reads
- 統(tǒng)計參考基因組上面各條染色體的成功比對reads數(shù)量
$grep -v '^@' tmp.sam | awk '!match($6,"*"){print($3);}' | sort | uniq -c
18995 gi|9626243|ref|NC_001416.1|
- 篩選出原始fq序列里面有N的比對情況
$grep -v '^@' tmp.sam | awk 'match($10,"N"){print}'
$grep -v '^@' tmp.sam | awk 'match($10,"N"){print}'| wc -l
12934
$grep -v '^@' tmp.sam | awk 'match($10,"N"){print}' | \
head -10 | cut -f 10 | grep "N"
- 篩選出原始fq序列里面有N毁腿，但是比對的時候卻是完全匹配的情況
$grep -v '^@' tmp.sam | awk 'match($6,length($10)){print}'
$grep -v '^@' tmp.sam | awk 'match($6,length($10)){print}'| wc -l
17095
$grep -v '^@' tmp.sam | awk 'match($6,length($10)){print}'| head -10 | less -SN
- sam文件里面的頭文件行數(shù)
$grep '^@' tmp.sam | wc -l
3
- sam文件里每一行的tags個數(shù)一樣嗎
$grep -v '^@' tmp.sam |awk '{print(NF-12+1)}'| sort | uniq -c
457 1
548 2
579 8
18416 9
- sam文件里每一行的tags個數(shù)分別是多少個
$grep -v '^@' tmp.sam |awk '{print(NF-12+1)}'
$grep -v '^@' tmp.sam |awk '{print(NF-12+1)}'|head -1
9
$grep -v '^@' tmp.sam |head -1 |\
awk '{for(i=12;i<NF;i++){printf("%s\n",$i)}}' | less -SN
013-2.png
- sam文件里記錄的參考基因組染色體長度分別是墅冷？
$grep "LN:" tmp.sam
@SQ SN:gi|9626243|ref|NC_001416.1| LN:48502
- 找到比對情況有insertion情況的
$grep -v '^@' tmp.sam | awk '{if (match($6,"[Ii]")) print;}'
$grep -v '^@' tmp.sam | awk '{if (match($6,"[Ii]")) print;}'| wc -l
66
$grep -v '^@' tmp.sam | awk '{if (match($6,"[Ii]")) print;}'|\
head -10 | less -SN
015.png
- 找到比對情況有deletion情況的
$grep -v '^@' tmp.sam | awk '{if (match($6,"[Dd]")) print;}'
$grep -v '^@' tmp.sam | awk '{if (match($6,"[Dd]")) print;}'| wc -l
1843
$grep -v '^@' tmp.sam | awk '{if (match($6,"[Dd]")) print;}'|\
head -10 | less -SN
17）取出位于參考基因組某區(qū)域的比對記錄，比如 5013到50130 區(qū)域
$grep -v '^@' tmp.sam | awk '{if ($4>=5013 && $4<=50130) print;}'
$grep -v '^@' tmp.sam | awk '{if ($4>=5013 && $4<=50130) print;}'| wc -l
17645
$grep -v '^@' tmp.sam | awk '{if ($4>=5013 && $4<=50130) print;}'|\
head -10 | less -SN
- 把sam文件按照染色體以及起始坐標(biāo)排序
$grep -v '^@' tmp.sam | sort -k3,4 -V
$grep -v '^@' tmp.sam | sort -k3,4 -V |head -20 |less -SN
$grep -v '^@' tmp.sam | sort -k3,4 -V |tail -20 |less -SN
head -20結(jié)果
tail - 20結(jié)果
- 找到 102M3D11M 的比對情況宛瞄，計算其reads片段長度。
$grep -v '^@' tmp.sam | awk '{if($6=="102M3D11M") print length($10);}'
113
- 安裝samtools軟件后使用samtools軟件的各個功能嘗試把上述題目重新做一遍交胚。
另開一篇文章寫作業(yè)份汗。
