2019-11-09 BAM和SAM格式文件shell练习

题目来源：http://www.bio-info-trainee.com/3578.html
首先使用bowtie2软件自带的测试数据生成sam/bam文件，代码如下（这里bowtie2是直接使用可执行文件，也可以source激活后直接使用）：

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

命令解释（参考生信技能树）：
bowtie的使用分为两步，首先建立索引，之后进行比对。
建索引：
基本命令：bowtie2-build your-fastq-file.fa your-index-name
example中的read已经建立好了索引，索引位置及前缀为../index/lambda_virus
比对基本格式：
单端测序：Bowtie –x tmp –U（表示单端） reads.fa –S hahahhha.sam
双端测序：Bowtie –x tmp -1 reads1.fa -2 reads2.fa –S hahahha.sam

../../bowtie2 还没有将bowtie2写入环境变量，因此使用bowtie2需要写明全路径
../index/lambda_virus 表示建立好的索引文件
-S 后接输出文件名
总结：将reads1.fa及reads2.fa 两条序列比对至索引文件../index/lambda_virus上比对后产生tmp.sam文件

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

less -SN tmp.sam | wc -l
20003

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

less -SN tmp.sam | sed '1,3d' | cut -f 2|sort -n | uniq -c
     24 65
    165 69
    153 73
    213 77
      2 81
   4650 83
    136 89
   4516 99
    125 101
     16 113
     24 129
    153 133
    165 137
    213 141
   4516 147
    125 153
      2 161
   4650 163
    136 165
     16 177

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

less -SN tmp.sam | grep -v '^@' |cut -f 1-10|grep -v 'M' |cut -f 10 |head
GCGTNGTACNNGNNGATGNTTACGACTAANACTTNNNNNTNCNGNNT
ANTGGNTNAGCGGNCGNAGNGNTCNGNAGCNANNAANCTGNTATGNTNNNTNNATCNNNCNNNNGTGNNAANNNCCAGNNNNGGCNNNNCATACANNNTNNNCNTNNNTACANNATANNCGNNATTNNACNGTNNCGNGNNCCGNGTNNNNN
NNTNNNNNAGNTNAANANANGTNTTGNGAANCTNNANACNNCTTTTNNNNTATGNTNNTNAGCCTAG
TTGNNCNNATATNNTAGATCGATATCTTTATNTCGAGCAAAGCAAA
NNNGANNNTGCTCGNCAANNNNANNNNTTTCCTTNTNCGGNCCNGNAATNATNACANNATGNNGNAANCGNNGNTTNTTNTNCAGGNCCTCTGNNGCNNANGCNCANGNNNTNTNNNNCNTNGNNGNCNNNCCGNAATNNNCTGNNGNCNTNGTNTCTGNCTNNNAGNNGNCGAANCTANNGANCCATANNTN
TTCGGNNCGCNNNNAAANAGGCTGAGCACGGTGTACGTCAGNCCGGAAAAGTGC
NNNGNTGTACCGGCTGTCTGGTATGTATGAGTTTGTGGTGAATAATGCCCCTGAACAGACAGAGGACGCCGGGCCCGCAGAGCCTGTTTCTGCGGGAAAGTGTTCGACGGTGAGCTGAGNTTTGCCCTGAAACTGG

less -SN tmp.sam |grep -v '^@' | cut -f 6 | grep -v 'M' | wc -l
1005

less -SN tmp.sam |grep -v '^@' |awk ' $6 =="*"  {print $0}' |wc -l
1005

less -SN tmp.sam |grep -v '^@' |awk '{if ( $6 =="*")  print $0}' |wc -l

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

less -SN tmp.sam |grep -v "^@" |awk '{if ($6=="*") print}' |cut -f1|sort -n |uniq -c|grep -w 2 #显示两端都没有比对上的id， -w参数表示精确匹配。

less -SN tmp.sam |grep -v "^@" |awk '{if ($6=="*") print}' |cut -f1|sort -n |uniq -c|grep -w 1

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

less -SN tmp.sam |grep -v "^@" |awk '{if ($5 > 30) print $0}' |cut -f5|sort -n |uniq -c
  1056 40
 17576 42

sam文件中，第一列为ID，第二列为比对类型，第三列是参考基因组信息，第四列为坐标，第五列为质量值，第6列为比对信息，第十列为碱基序列

筛选出比对成功，但是并不是完全匹配的序列

cat tmp.sam |grep -v "^@" |cut -f 6,10|grep [IDPSH]|wc -l

M：alignment match (can be a sequence match or mismatch)

表示read可mapping到第三列的序列上，则read的碱基序列与第三列的序列碱基相同，表示正常的mapping结果，M表示完全匹配，但是无论reads与序列的正确匹配或是错误匹配该位置都显示为M

I：insertion to the reference

表示read的碱基序列相对于第三列的RNAME序列，有碱基的插入

D：deletion from the reference

表示read的碱基序列相对于第三列的RNAME序列，有碱基的删除

N：skipped region from the reference

表示可变剪接位置

P：padding (silent deletion from padded reference)

S：soft clipping (clipped sequences present in SEQ)

H：hard clipping (clipped sequences NOT present in SEQ)

clipped均表示一条read的序列被分开，之所以被分开，是因为read的一部分序列能匹配到第三列的RNAME序列上，而被分开的那部分不能匹配到RNAME序列上。

"="表示正确匹配到序列上

"X"表示错误匹配到序列上

而H只出现在一条read的前端或末端，但不会出现在中间，S一般会和H成对出现，当有H出现时，一定会有一个与之对应的S出现

例如：

162M89S

162H89M

149M102S

149H102M

40S211M

20M1D20M211H

S可以单独出现，而H必须有与之对应的S出现时才可能出现，不可在相同第一列的情况下单独出现

N：如果是mRNA-to-genome，N出现的位置代表内含子，其它比对形式出现N时则没有具体解释

M/I/S/=/X：这些数值的加和等于第10列SEQ的长度

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

less -SN tmp.sam | grep -v '^@' | awk '{if( $9 > 1250 || $9 < -1250) print $0}' |wc -l
80

less -SN tmp.sam | grep -v '^@' | awk ' $9 > 1250 || $9 < -1250 { print $0}' |wc -l
80

统计参考基因组上面各条染色体的成功比对reads数量.


less -SN tmp.sam | grep -v '^@' | cut -f 3 | sort |uniq -c
      426 *
  19574 gi|9626243|ref|NC_001416.1|

筛选出原始fq序列里面有N的比对情况

less -SN tmp.sam | grep -v '^@' |cut  -f  10 | grep 'N' |wc -l
12934

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

less -SN tmp.sam | grep -v '^@' |cut  -f  6,10 |grep 'N'|cut -f 1|grep -v [IDNPSH] |grep 'M'|wc -l
10632

less -SN tmp.sam | grep -v '^@' |awk ' $10 ~ "N" {print $0}' | awk '$6 !~  "[IDNPSH]" {print $0}'| awk '$6 !~ "*" {print $0}' |wc -l
10632

less -SN tmp.sam | grep -v '^@' |awk ' $10 ~ "N" {print $0}' | awk '$6 !~  "[IDNPSH]" {print $0}'| awk '$6 != "*" {print $0}' |wc -l
10632

less -SN tmp.sam | grep -v  '^@' | awk '{if ($10 ~ "N") print $0}' | awk '{ if ($6 !~ "[IDNPSH]")  print $0} '| awk '{if  ($6 != "*") print $0}'  | wc -l
10632

less -SN tmp.sam | grep -v  '^@' | awk '{if ($10 ~ "N") print $0}' | awk '{ if ($6 !~ "[IDNPSH]")  print $0} '| awk '{if  ($6 != "\*") print $0}'  | wc -l
awk: cmd. line:1: warning: escape sequence `\*' treated as plain `*'  # awk命令提示"*"字符不需要转义符
10632

less -SN tmp.sam | grep -v '^@' | awk '$10 ~ "N" && $6 !~ "[IDNPSH]"  && $6 != "*" {print $0}'|wc -l
10632

less -SN tmp.sam | grep -v '^@'| awk '{if ($10 ~ "N" && $6 !~ "[IDNPSH]"  && $6 != "*") print $0}' |wc -l
10632

sam文件里面的头文件行数

less -SN tmp.sam | grep  '^@' |wc -l
3

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

less -SN tmp.sam | grep -v '^@' |cut -f 12- | less -SN

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

less -SN tmp.sam | grep -v '^@' |cut -f 12| wc -l
20000
less -SN tmp.sam | grep -v '^@' |cut -f 20| wc -l
20000 #第12列和行数与第20列一致，说明每一行的tags数量相同，均为9个

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

less -SN tmp.sam |grep '^@' |grep 'LN'
@SQ SN:gi|9626243|ref|NC_001416.1|  LN:48502

找到比对情况有insertion情况的

less -SN tmp.sam | grep -v '^@' | awk ' $6 ~ "I"  {print $0}' | wc -l
66

找到比对情况有deletion情况的

less -SN tmp.sam | grep -v '^@' | awk ' $6 ~ "I"  {print $0}' | wc -l
1843

取出位于参考基因组某区域的比对记录，比如 5013到50130 区域

less -SN tmp.sam| grep -v '^@' | awk '$4 > 5013 && $4 < 50130 && $8 > 5013 && \
  $8 < 50130 && $6 != "*" { print $0}' |wc -l
17072

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

less -SN tmp.sam| grep -v '^@' | sort   -k 4 -n | awk '$4 != 0 {print $0}'|wc -l
19574

less -SN tmp.sam| grep -v '^@' | sort  -t ' ' -k 4 -n | awk '$4 != 0 {print $0}'|wc -l
19574

找到 102M3D11M 的比对情况，计算其reads片段长度

less -SN tmp.sam| grep -v '^@' | grep -w  '102M3D11M'| awk '{print length($10)}'
113

安装samtools软件后使用samtools软件的各个功能尝试把上述题目重新做一遍。
其它概念题
高通量测序数据分析中，序列与参考序列的比对后得到的标准文件为什么文件？
sam文件是一种比对后的文件，能直接查看吗？
Sam/Bam文件分为两部分，一部分以@开头的是什么序列，另一部分是什么序列？
Sam文件除头文件，以什么符号分割文本的，比对部分信息一行是多少列？你能用awk计算列数吗？
Sam/Bam文件的@SQ开头的行是什么？你能生成一个文本，两列，一列是参考基因组染色体/sca id，一列是长度吗？
Sam文件的比对位置是从1还是0开始的？
常见CIGAR 字符串各字母代表的意义
比对质量的数字是哪一列？越大比对质量越好还是越小越好？
Sam文件的flag是第几列，数字代表什么意义？是怎么计算来的？
Sam文件怎么转bam文件？用什么指令？为什么要转换？
Bam文件查看用什么指令？如果需要查看头文件需要增加参数？
Bam文件为什么要排序？排序后的bam和未排序的bam头文件和chr position列有什么区别？
Bam文件建索引的指令是什么？
Bam文件可视化用什么工具？查看时需要建立索引吗？
Bam文件统计reads比对情况用samtools的哪一个子命令？
SE测序和PE测序的所比对后得到的sam文件的区别在哪里？
Bam能用gzip再压缩吗？
Sam文件通常由哪些比对软件得到的
Sam/Bam文件能转成fastq文件吗？
有时候不能通过文件名的后缀来区别是sam还是bam，你是怎么区别的？