生信学习——sam和bam格式文件的shell小练习（附详细答案解读）

sam的全称是：sequence alignment/map format。bam是sam的二进制文件，sam文件包括标头注释部分和比对结果部分。

本文主要介绍如何使用Linux中的shell命令来处理这种文件。

视频教程：https://www.bilibili.com/video/BV1ds411g7eg?p=13
题目原文：http://www.bio-info-trainee.com/3578.html
格式说明：http://samtools.github.io/hts-specs/

准备练习文件

首先使用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
# 生成sam文件
../../bowtie2 -x ../index/lambda_virus -1 reads_1.fq -2 reads_2.fq > tmp.sam

# 使用samtools将sam文件转成bam文件
# samtools view -bS tmp.sam >tmp.bam

tmp.sam文件内容：

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

# 思路：去除开头前三行，利用第一列去重，再统计
cat tmp.sam | grep -v '^@' | cut -f 1 | uniq | wc -l

###
# 输出结果
10000

注：双端测序(Paired-end和Mate-pair)，是对一个长的序列测得其两端的序列。两端的序列形成"一对"，中间的距离叫插入长度（insert length）。

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

# -n是按照数字大小排序，-r是以相反顺序，-k是指定需要排序的列数
cat tmp.sam | grep -v '^@' | cut -f 2 | sort | uniq -c | sort -nrk 1

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

Decoding SAM flags：http://broadinstitute.github.io/picard/explain-flags.html
example：

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

注：第六列为简要比对信息表达式，如果值为“*”则表示比对失败。

# 答案1
cat tmp.sam | grep -v '^@' | cut -f 6 | grep '*' | wc
# 答案2
cat tmp.sam | grep -v '^@' | awk '{if($6=="*") print}'| wc

序列特征：含有较多的N碱基

​

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

# 单端失败 ID在比对失败的reads列表中只出现一次
cat tmp.sam | grep -v '^@' | awk '{if($6=="*") print$1}' | sort | uniq -c | grep -w 1

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

​

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

# print$0与print 都是输出全部信息
cat tmp.sam | grep -v '^@' | awk '{if($5 > 30) print$0}' | less -S

​

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

# 第六列出现IDNSHPX其中一个，则表示未完全匹配
# ~ 匹配正则表达式
cat tmp.sam | grep -v '^@' | awk '{if($6 ~ "[IDNSHPX]") print}' | less -S

​

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

cat tmp.sam | grep -v '^@' | awk '{if($9>1250 || $9<-1250) print}' | wc

###
#输出结果
     80    1600   25902

​

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

cat tmp.sam | grep -v '^@' | cut -f 3 | sort | uniq -c

###
# 输出结果
    426 *
  19574 gi|9626243|ref|NC_001416.1|

成功比对数量为：19574
​

9. 筛选出原始fq序列里面有N的比对情况（看第10列）

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

###
#输出结果
  12934  251331 4881351

​

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

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

###
#输出结果
  10632  212320 3994276

​

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

grep -c '^@' tmp.sam 

###
#输出结果
3

​

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

# tag在11列之后
grep -v '^@' tmp.sam | cut -f 12- tmp.sam | less -S

从图片中可以看到每行的tags个数不是全都一样的。
​

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

#  NF 表示字段数，在执行过程中对应于当前的字段数。
grep -v '^@' tmp.sam | cut -f 12- | awk '{print NF}' | sort | uniq -c

###
#输出结果
    457 1
    548 2
    579 8
  18416 9

​

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

# 头部信息中的LN即为染色体长度
head -3 tmp.sam 

​

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

# insertion代表第六列包含“I”
grep -v '^@' tmp.sam | awk '{if($6 ~ "I") print}' | less -S

​

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

# deletion代表第六列包含“D”
grep -v '^@' tmp.sam | awk '{if($6 ~ "D") print}' | less -S

​

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

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

​

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

# if($4 != 0) 去除无法比对的染色体
# -nk 1.2 从第一个域的第二个字符开始排序
# -nk 4 根据第四列进行排序
grep -v '^@' tmp.sam | awk '{if($4 != "0") print}' | sort -nk 1.2 -nk 4 | less -SN

​

19. 找到 102M3D11M 的比对情况，计算其reads片段长度。（看第10列）

grep '102M3D11M' tmp.sam | awk '{print length($10)}'

###
# 输出结果
113

​

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

先挖个坑，等后期学习samtools的使用时，过来填坑。
​

其它概念题

优质答案：https://hc1023.github.io/2019/09/07/linux-sam-bam/

1. 高通量测序数据分析中，序列与参考序列的比对后得到的标准文件为什么文件？

为sam文件

2. sam文件是一种比对后的文件，能直接查看吗？

可以使用cat head less等命令直接查看

3. Sam/Bam文件分为两部分，一部分以@开头的是什么序列，另一部分是什么序列？

以@开头的是标头注释序列，另一部分是比对结果序列

4. Sam文件除头文件，以什么符号分割文本的，比对部分信息一行是多少列？你能用awk计算列数吗？

11个必须字段和1个可选字段，字段之间用tag分割

# 第13题
grep -v '^@' tmp.sam | cut -f 12- | awk '{print NF}' | sort | uniq -c

5. Sam/Bam文件的@SQ开头的行是什么？你能生成一个文本，两列，一列是参考基因组染色体/sca id，一列是长度吗？

参考序列说明

grep '^@SQ' tmp.sam | cut -f 2,3
# 输出结果
SN:gi|9626243|ref|NC_001416.1|	LN:48502

6. Sam文件的比对位置是从1还是0开始的？

从1开始的

7. 常见CIGAR 字符串各字母代表的意义

8. 比对质量的数字是哪一列？越大比对质量越好还是越小越好？

第5列，越大说明该read比对到参考基因组上的位置越准确

9. Sam文件的flag是第几列，数字代表什么意义？是怎么计算来的？

第2列，数字意义如图，如何计算参考第2题。

10. Sam文件怎么转bam文件？用什么指令？为什么要转换？

samtools view -bS tmp.sam >tmp.bam

bam为二进制文件，占用内存小

11. Bam文件查看用什么指令？如果需要查看头文件需要增加参数？

samtools view tmp.bam | less -S
samtools view -H tmp.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，你是怎么区别的？

注：由于本人对samtools、bam文件还不熟悉，所以后面几题请直接查看大佬的优质答案。