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

原文链接

首先使用bowtie2软件自带的测试数据生成sam/bam文件，代码如下：

mkdir -p ~/biosoft
cd ~/biosoft
wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.3.4.3/bowtie2-2.3.4.3-linux-x86_64.zip 
unzip bowtie2-2.3.4.3-linux-x86_64.zip 
cd ~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads
../../bowtie2 -x ../index/lambda_virus -1 reads_1.fq -2 reads_2.fq > tmp.sam
# samtools view -bS tmp.sam >tmp.bam
#利用以上命令生成：tmp.sam 与tmp.bam两个文件

LINUX练习题

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

cat tmp.sam | grep -v '^@'|cut -f1|sort |uniq -c|wc -l #去掉三行以@开头的（头文件）
10000

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

 cat tmp.sam |grep -v "^@"|awk '{print $2}'|sort -n|uniq -c|less -N
1      24 65
      2     165 69
      3     153 73
      4     213 77
      5       2 81
      6    4650 83
      7     136 89
      8    4516 99
      9     125 101
     10      16 113
     11      24 129
     12     153 133
     13     165 137
     14     213 141
     15    4516 147
     16     125 153
     17       2 161
     18    4650 163
     19     136 165
     20      16 177
#or
cat tmp.sam |grep -v "^@"|cut -f 2|sort |uniq -c|less -N

cat tmp.sam | grep -v '^@' | cut -f2|sort|uniq -c|sort -k1,1nr|less -SN

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

cat tmp.sam |grep -v '^@'|awk '{if($6 =="*") print}'|wc #第6列为" * " 表示比对失败
   1005   12608  255140

cat tmp.sam |grep -v '^@'|awk '{if($6 =="*") print$10}' # 第10列为序列，序列含有较多的“N”碱基

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

#单端失败
cat tmp.sam |grep -v '^@'|awk '{if($6 =="*") print$1}'|sort -n |uniq -c |grep -w 1 # 序列ID为“$1”，grep -w 1表示ID出现一次

#双端失败
cat tmp.sam |grep -v '^@'|awk '{if($6 =="*") print$1}'|sort -n |uniq -c |grep -w 2

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

cat tmp.sam |grep -v '^@'|awk '{if($5>30) print}' |wc #统计出
  18632  372088 6661992
cat tmp.sam |grep -v '^@'|awk '{if($5>30) print}' |less -SN > 大于Q30.txt

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

cat tmp.sam |grep -v '^@'|awk '{if($6!="*") print$6}'|grep [IDNXSHP]
“M”表示 match或 mismatch；
“I”表示 insert
“D”表示 deletion
“N”表示 skipped（跳过这段区域）
“S”表示 soft clipping（被剪切的序列存在于序列中）
“H”表示 hard clipping（被剪切的序列不存在于序列中）
“P”表示 padding
“=”表示 match
“X”表示 mismatch（错配，位置是一一对应的）

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

cat tmp.sam |grep -v '^@'|awk '{if($9>1250) print}'|less -SN #视频里面$7是错的

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

cat tmp.sam |grep -v '^@'|cut -f 3 |sort -u
*
gi|9626243|ref|NC_001416.1|

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

cat tmp.sam |grep -v '^@'|awk '{if($10 ~"N")print}'|less -SN #匹配

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

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

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

cat tmp.sam |grep '^@'|wc
      3      19     269

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

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

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

代码块

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

cat tmp.sam | grep '^@' | grep 'LN'

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

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

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

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

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

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

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

cat tmp.sam | grep -v '^@'|sort -k4|less -SN

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

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

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

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

后续更新

其它概念题

1. 高通量测序数据分析中，序列与参考序列的比对后得到的标准文件为什么文件？
2. sam文件是一种比对后的文件，能直接查看吗？
3. Sam/Bam文件分为两部分，一部分以@开头的是什么序列，另一部分是什么序列？
4. Sam文件除头文件，以什么符号分割文本的，比对部分信息一行是多少列？你能用awk计算列数吗？
5. Sam/Bam文件的@SQ开头的行是什么？你能生成一个文本，两列，一列是参考基因组染色体/sca id，一列是长度吗？
6. Sam文件的比对位置是从1还是0开始的？
7. 常见CIGAR 字符串各字母代表的意义
8. 比对质量的数字是哪一列？越大比对质量越好还是越小越好？
9. Sam文件的flag是第几列，数字代表什么意义？是怎么计算来的？
10. Sam文件怎么转bam文件？用什么指令？为什么要转换？
11. Bam文件查看用什么指令？如果需要查看头文件需要增加参数？
12. Bam文件为什么要排序？排序后的bam和未排序的bam头文件和chr position列有什么区别？
13. Bam文件建索引的指令是什么？
14. Bam文件可视化用什么工具？查看时需要建立索引吗？
15. Bam文件统计reads比对情况用samtools的哪一个子命令？
16. SE测序和PE测序的所比对后得到的sam文件的区别在哪里？
17. Bam能用gzip再压缩吗？
18. Sam文件通常由哪些比对软件得到的
19. Sam/Bam文件能转成fastq文件吗？
20. 有时候不能通过文件名的后缀来区别是sam还是bam，你是怎么区别的？