生信格式SAM、BAM

资料推荐

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
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

学徒培养：https://mp.weixin.qq.com/s/3jw3_PgZXYd7FomxEMxFmw

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
【社区内容提示】社区部分内容疑似由AI辅助生成，浏览时请结合常识与多方信息审慎甄别。
平台声明：文章内容（如有图片或视频亦包括在内）由作者上传并发布，文章内容仅代表作者本人观点，简书系信息发布平台，仅提供信息存储服务。

相关阅读更多精彩内容

  • 八爪鱼采集器
    八爪鱼采集器 2018-02-08星夜读书区 今天，星夜小编为大家讲解八爪鱼采集器是什么、八爪鱼采集器怎么用和八爪...
    最美有品阅读 2,804评论 0赞 2
  • 2016.07.10
    1、情绪管理：相处愉快，没有情绪冲突 2、亲子陪伴:今天全家去看望刚生病出院的太奶奶。我帮爷爷做饭打下手，小朋友自...
    Emmaluo阅读 217评论 0赞 0
  • 谁之罪
    今天去看了我们学校的一场话剧《谁之罪》，讲的是一个刚刚获得奥斯卡最佳导演的导演布鲁斯，热衷在电影中宣言暴力。谁知，...
    萌小诗zzz阅读 479评论 0赞 0

友情链接更多精彩内容