资料推荐
- 生信菜鸟团的浅谈SAM格式
- 生信媛的高通量数据分析必须知道的山姆大叔(SAM)
- 生信技能树系列视频数据格式,以及SAM、BAM练习题讲解
- 生信菜鸟团的sam和bam格式文件的shell小练习
sam和bam格式文件的shell小练习
1) 统计共多少条reads(pair-end reads这里算一条)参与了比对参考基因组
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat tmp.sam | grep -v '^@'|wc
20000 391929 7049181
# 左右两个序列算作一条，所以为10000
2) 统计共有多少种比对的类型(即第二列数值有多少种)及其分布。
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat tmp.sam | grep -v '^@' | cut -f2|sort|uniq -c|sort -k1,1nr
4650 163
4650 83
4516 147
4516 99
213 141
213 77
165 137
165 69
153 133
153 73
136 165
136 89
125 101
125 153
24 129
24 65
16 113
16 177
2 161
2 81
- picard sam flag:https://broadinstitute.github.io/picard/explain-flags.html
picard sam flag
3）筛选出比对失败的reads，看看序列特征。
# 第6列的* 代表为比对失败
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat tmp.sam | grep -v '^@' |awk '{if($6=="*")print}'|wc
1005 12608 255140
4) 比对失败的reads区分成单端失败和双端失败情况，并且拿到序列ID
# 单端失败
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat tmp.sam | grep -v '^@' |awk '{if($6=="*")print $1}'|sort|uniq -c|grep -w 1
# 双端失败
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat tmp.sam | grep -v '^@' |awk '{if($6=="*")print $1}'|sort|uniq -c|grep -w 2
5) 筛选出比对质量值大于30的情况（看第5列）
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat tmp.sam | grep -v '^@' |awk '{if($5>30)print}'|wc
18632 372088 6662664
6) 筛选出比对成功，但是并不是完全匹配的序列
“M”表示 match或 mismatch；
“I”表示 insert
“D”表示 deletion
“N”表示 skipped（跳过这段区域）
“S”表示 soft clipping（被剪切的序列存在于序列中）
“H”表示 hard clipping（被剪切的序列不存在于序列中）
“P”表示 padding
“=”表示 match
“X”表示 mismatch（错配，位置是一一对应的）
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat tmp.sam | grep -v '^@' |awk '{if($6!="*")print$6}'|grep "[IDNSHPX]"|wc
1900 1900 18522
7) 筛选出inset size长度大于1250bp的 pair-end reads
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat tmp.sam | grep -v '^@' |awk '{if($7>1250)print}'|less -S
8) 统计参考基因组上面各条染色体的成功比对reads数量
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat tmp.sam | grep -v '^@' |cut -f3|sort -u
*
gi|9626243|ref|NC_001416.1|
9) 筛选出原始fq序列里面有N的比对情况
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat tmp.sam | grep -v '^@'|awk '{if($10~N)print}'|wc
20000 391929 7049181
10) 筛选出原始fq序列里面有N，但是比对的时候却是完全匹配的情况
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat tmp.sam | grep -v '^@'|awk '{if($10 ~ N)print}'|awk '{if($6 !~ "[IDNSHP]")print}'|awk '{if($6!="*")print}'|less -S
11) sam文件里面的头文件行数
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ grep '^@' tmp.sam|wc
3 19 262
12) sam文件里每一行的tags个数一样吗;13) sam文件里每一行的tags个数分别是多少个
14) sam文件里记录的参考基因组染色体长度分别是？
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ head tmp.sam | grep 'LN'
@SQ SN:gi|9626243|ref|NC_001416.1| LN:48502
15) 找到比对情况有insertion情况的
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat tmp.sam | grep -v '^@'|awk '{if($6~I)print}'|less -S
16) 找到比对情况有deletion情况的
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat tmp.sam | grep -v '^@'|awk '{if($6~D)print}'|less -S
17）取出位于参考基因组某区域的比对记录，比如 5013到50130 区域
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat tmp.sam | grep -v '^@'|awk '{if($4>5013 && $4 <50130)print}'|less -S
18) 把sam文件按照染色体以及起始坐标排序
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ cat tmp.sam | grep -v '^@'|awk '{print $4}'|sort -n
# 还有点问题
19) 找到 102M3D11M 的比对情况，计算其reads片段长度
vip39@VM-0-15-ubuntu:~/test/bowtie2-2.3.4.3-linux-x86_64/example/reads$ grep 102M3D11M
tmp.sam |cut -f 10|wc
1 1 114
# 所以就114咯
20) 安装samtools软件后使用samtools软件的各个功能尝试把上述题目重新做一遍。
vip39@VM-0-15-ubuntu:~$ source miniconda3/bin/activate
(base) vip39@VM-0-15-ubuntu:~$ cd src/
(base) vip39@VM-0-15-ubuntu:~/src$ conda install samtools=1.7 y
生信技能树公益视频合辑：学习顺序是linux，r，软件安装，geo，小技巧，ngs组学！
请猛戳下面链接：
B站链接：https://m.bilibili.com/space/338686099
YouTube链接：https://m.youtube.com/channel/UC67sImqK7V8tSWHMG8azIVA/playlists
生信工程师入门最佳指南：https://mp.weixin.qq.com/s/vaX4ttaLIa19MefD86WfUA
