【生信技能树】sam和bam格式文件的shell小练习

LINUX练习题

1) 统计共多少条reads(pair-end reads这里算一条)参与了比对参考基因组

$grep -v  '^@' tmp.sam |cut -f 1|sort -V | uniq -c | wc -l
10000 

2. 统计共有多少种比对的类型(即第二列数值有多少种)及其分布。

$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

4. 比对失败的reads区分成单端失败和双端失败情况，并且拿到序列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

5. 筛选出比对质量值大于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

6. 筛选出比对成功，但是并不是完全匹配的序列

7. 筛选出inset size长度大于1250bp的 pair-end reads

8. 统计参考基因组上面各条染色体的成功比对reads数量

$grep -v '^@' tmp.sam | awk '!match($6,"*"){print($3);}' | sort | uniq -c
  18995 gi|9626243|ref|NC_001416.1|

9. 筛选出原始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"

10. 筛选出原始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

11. sam文件里面的头文件行数

$grep '^@' tmp.sam | wc -l
3

12. sam文件里每一行的tags个数一样吗

$grep -v  '^@' tmp.sam |awk '{print(NF-12+1)}'| sort | uniq -c
    457 1
    548 2
    579 8
  18416 9

13. sam文件里每一行的tags个数分别是多少个

$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

14. sam文件里记录的参考基因组染色体长度分别是？

$grep "LN:" tmp.sam
@SQ     SN:gi|9626243|ref|NC_001416.1|  LN:48502

15. 找到比对情况有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

16. 找到比对情况有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）取出位于参考基因组某区域的比对记录，比如 5013到50130 区域

$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

18. 把sam文件按照染色体以及起始坐标排序

$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结果
tail - 20结果

19. 找到 102M3D11M 的比对情况，计算其reads片段长度。

$grep -v  '^@' tmp.sam | awk '{if($6=="102M3D11M") print length($10);}'
113

20. 安装samtools软件后使用samtools软件的各个功能尝试把上述题目重新做一遍。

另开一篇文章写作业。