生信学习——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