sam和bam格式文件的shell小练习

经过一段时间的学习，不断的逛生信菜鸟团和生信技能树对RNA-SEQ有个基本的了解，于是静下心来做下作业。

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

第一步先按照要求下载数据

less -SN tmp.sam|grep -v '^@'|wc

结果显示为20000，但由于是双端测序，所以答案为10000.

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

less -SN tmp.sam|grep -v '^@'|cut -f 2|sort |uniq -c

125 101
     16 113
     24 129
    153 133
    165 137
    213 141
   4516 147
    125 153
      2 161
   4650 163
    136 165
     16 177
     24 65
    165 69
    153 73
    213 77
      2 81
   4650 83
    136 89
   4516 99

3.筛选出比对失败的reads，看看序列特征。

第6列的* 代表为比对失败

less -SN tmp.sam|grep -v '^@'|awk '{if ($6=="*")print}'|wc

结果为1005

4. 比对失败的reads区分成单端失败和双端失败情况，并且拿到序列ID

序列ID为第一列

less -SN tmp.sam|grep -v '^@'|awk '{if ($6=="*")print $1}'|sort|unqi -c|grep -w 1
为一端没有比对上。同理
less -SN tmp.sam|grep -v '^@'|awk '{if ($6=="*")print $1}'|sort|unqi -c|grep -w 2
为二端没有比对上。

5.筛选出比对质量值大于30的情况（看第5列）

less -SN tmp.sam|grep -v '^@'|awk '{if ($5>30)print }'|wc

结果显示为18632.

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

完全匹配为M，比对成功意味着第6列不是*，比对成功，但不是完全匹配，有IDNSHPX这几种情况。

less -SN tmp.sam|grep -v '^@'|awk '{if ($6!="*")print $6}'| grep "[IDNSHPX]"|wc

结果显示为1900.

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

less -SN tmp.sam|grep -v '^@' |awk '{if($7>1250)print}'|less -S

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

cut -f 3 tmp.sam|sort

9.筛选出原始fq序列里面有N的比对情况

less -SN tmp.sam|grep -v '^@' |awk '{if($10~N)print}'|less -S|wc

10.筛选出原始fq序列里面有N，但是比对的时候却是完全匹配的情况

less -SN tmp.sam|grep -v '^@' |awk '{if($10~N)print}'|awk '{if($6!~"[IDNSHP]")print}'|awk '{if($6!~"*") print}'|less -SN

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

less -SN tmp.sam|grep - '^@'|wc

结果显示为3

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

cat tmp.sam | grep -v '^@' |cut -f 12- |less -S

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

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

14题和13题是一块的，其长度为48052.

15.找到比对情况有insertion情况的

less  tmp.sam | grep -v '^@'|awk '{if($6~I)print}'|less -S

16.找到比对情况有deletion情况的

思路是一样的，

less  tmp.sam | grep -v '^@'|awk '{if($6~D)print}'|less -S

17.取出位于参考基因组某区域的比对记录，比如 5013到50130 区域

less  tmp.sam | grep -v '^@'|awk '{if($4>5013 && $4 <50130)print}'|less -S

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

less  tmp.sam | grep -v '^@'|awk '{print $4}'|sort -n

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

grep  102M3D11M tmp.sam |cut -f 10|wc