先使用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两个文件
1.统计共多少条reads(pair-end reads这里算一条)参与了比对参考基因组
cat tmp.sam | grep -v '^@'|cut -f1|sort |uniq -c|wc -l #去掉三行以@开头的（头文件）
10000
2.统计共有多少种比对的类型(即第二列数值有多少种)及其分布。
cat tmp.sam | grep -v '^@'|cut -f2|sort |uniq -c|less -N
#或
cat tmp.sam | grep -v '^@'|awk '{print $2}'|sort |uniq -c|less -N
>
1 125 101
2 16 113
3 24 129
4 153 133
5 165 137
6 213 141
7 4516 147
8 125 153
9 2 161
10 4650 163
11 136 165
12 16 177
13 24 65
14 165 69
15 153 73
16 213 77
17 2 81
18 4650 83
19 136 89
20 4516 99
3.筛选出比对失败的reads，看看序列特征。
cat tmp.sam |grep -v '^@'|awk '{if($6 =="*") print}'|less -S #第6列为" * " 表示比对失败
1005 12608 255140
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}' |less -SN
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
8.统计参考基因组上面各条染色体的成功比对reads数量
cat tmp.sam |grep -v '^@'|cut -f 3 |sort -u
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 -l
12.sam文件里每一行的tags个数一样吗
cat tmp.sam | grep -v '^@' |cut -f 12- |less -S
13.sam文件里每一行的tags个数分别是多少个
###这里有 问题，不是用空格分开的，所以其实答案是错的
grep -v '^@' tmp.sam |cut -f 12-|awk '{split($0,a, " ");print length(a)}' |sort|uniq -c
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 -m
网友评论