用上一节下载的bowtie2软件中自带的测试数据生成sam和bam文件先：
# 下载bowtie2软件先
wget -c https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.3.5.1/bowtie2-2.3.5.1-linux-x86_64.zip
unzip bowtie2-2.3.5.1-linux-x86_64.zip
# 找到示例数据
cd bowtie2-2.3.5.1-linux-x86_64/example/reads
# 先生成sam文件，使用bowtie2命令，绝对路径调用，参数什么的DIY吧
/trainee/Jude/practice/bowtie2-2.3.5.1-linux-x86_64/bowtie2 -x /trainee/Jude/practice/bowtie2-2.3.5.1-linux-x86_64/example/index/lambda_virus -1 reads_1.fq -2 reads_2.fq > test.sam
# sam转换成bam
samtools view -bS test.sam > test.bam
一、统计共多少条reads(pair-end reads这里算一条)参与了比对参考基因组
# 这道题的命令很简单，但是在写命令之前花了好长的时间学习sam文件每一列的含义.下面命令的值是20000，因为是双端，所以有10000条pair-end reads参与了比对。
less -SN test.sam | grep -v "^@" | wc
二、统计共有多少种比对的类型(即第二列数值有多少种)及其分布
# 就是统计有多少种flag, 取出第二列，排序，然后合并计数
less -SN test.sam | grep -v "^@" | awk '{print $2}' | sort | uniq -c
三、筛选出比对失败的reads，看看序列特征
# 先搞明白什么叫作比对失败的reads：第6列 ( CIGAR )为*则表示比对失败,先计数
less -SN test.sam | grep -v "^@" | awk '{if($6=="*")print}'| wc
# 再看序列特征, 序列信息在第10例
less -SN test.sam | grep -v "^@" | awk '{if($6=="*")print}'| awk '{print$10}' | less -SN
其实这里有个疑问，按说第三列如果为*也表示比对失败（没有比对上），为什么最终的数值不一样呢···
四、比对失败的reads区分成单端失败和双端失败情况，并且拿到序列ID
# 考的是概念不是技术啊啊啊
less test.sam |grep -v '^@'|awk '{if($6 =="*") print$1}'|sort -n |uniq -c |grep -w 1
less test.sam |grep -v '^@'|awk '{if($6 =="*") print$1}'|sort -n |uniq -c |grep -w 2
五、筛选出比对质量值大于30的情况（看第5列）
# 送分题 : 第五列表示比对质量
less test.sam |grep -v '^@'|awk '{if($5 > 30) print}' | less -SN
less test.sam |grep -v '^@'|awk '{if($5 > 30) print}' | wc
六、筛选出比对成功，但是并不是完全匹配的序列
# 第一反应还是看第6列
less test.sam |grep -v '^@'|awk '{if($6!="*") print$6}'| grep [IDNXSHP]
# 纠结了很久，因为总觉得最后的展示结果应该是$10列，毕竟是要看序列嘛。后来一想，单纯看序列看不出来这些突变呃。。。
七、筛选出insert size长度大于1250bp的 pair-end reads
less test.sam | grep -v '^@'| awk '{if($9>1250) print}'| head | less -SN
八、统计参考基因组上面各条染色体的成功比对reads数量
# 成功比对后第3列会显示对应的参考序列信息，无法比对则显示*，所以取出第3列，排序然后计数
less test.sam | grep -v '^@'| cut -f 3 | sort -u | uniq -c
九、筛选出原始fq序列里面有N的比对情况
# 比对情况看第6列，不过首先需要筛选出有N的行
less test.sam | grep -v '^@'| awk '{if($10~N) print$6}' | less
十、筛选出原始fq序列里面有N，但是比对的时候却是完全匹配的情况
# 同上题思路
less test.sam | grep -v '^@'| awk '{if($10~N) print}' | awk '{if($6!~"[IDNXSHP*]") print}' | less -SN
发现卡壳主要是因为不熟悉awk的语法和SAM数据结构，这样的时候千万不要跟自己较劲，去看别人的答案和教程吧，因为不是逻辑推理的问题，而更像是公式，参考别人的，学会了就好。
十一、sam文件里面的头文件行数
# 头文件结构是以@符号打头
less test.sam | grep '^@' | wc
十二、sam文件里每一行的tags个数一样吗
# 搞清楚tags是啥，在什么地方先:在第12列及后面的列，为可选字段（optional fields)，格式如：TAG:TYPE:VALUE，其中TAG有两个大写字母组成，每个TAG代表一类信息，每一行一个TAG只能出现一次，TYPE表示TAG对应值的类型，可以是字符串、整数、字节、数组等。
less test.sam | grep -v '^@' | cut -f 12- | less -SN
# 肉眼可见的不一样
十三、sam文件里每一行的tags个数分别是多少个
# 跟上题一样, 由输出结果可知绝大部分为9个，也有2个，8个的
less test.sam | grep -v '^@' | cut -f 12- | awk -F '\t' '{print "columns:" NF }'
这道题做出来很开心，因为目前能搜到的其他人写的答案这道题都是空着的，哈哈哈~~
十四、sam文件里记录的参考基因组染色体长度分别是？
# 参考基因组染色体的长度信息记录在哪里？
less test.sam | grep '^@' | grep LN
十五、找到比对情况有insertion情况的
# 有insertion表示第6列里面有I
less test.sam | grep -v '^@'| cut -f 6 | grep I
十六、找到比对情况有deletion情况的
# 有deletion表示第6列里面有D
less test.sam | grep -v '^@'| cut -f 6 | grep D
十七、取出位于参考基因组某区域的比对记录，比如 5013到50130 区域
# 比对区域看第4列，表示比对的起始位置
less test.sam | grep -v '^@'| awk '{if($4>5013 && $4<50130) print}' | less -SN
十八、把sam文件按照染色体以及起始坐标排序
# 就看sort怎么用，首先以染色体排序：
less test.sam | grep -v '^@' | sort -k 3 | less -S
# 以其实坐标排序，坐标是数字，需-n参数
less test.sam | grep -v '^@' | sort -n -k 4 | less -S
查看结果可知如果按照升序排列的话最前面的都是星号（*），可以改为降序。
十九、找到 102M3D11M 的比对情况，计算其reads片段长度
# grep一下，然后cut或者awk一下
less test.sam | grep "102M3D11M" | cut -f 10 | wc
二十、安装samtools软件后使用samtools软件的各个功能尝试把上述题目重新做一遍
其它概念题
- 高通量测序数据分析中，序列与参考序列的比对后得到的标准文件为什么文件？
SAM文件
- sam文件是一种比对后的文件，能直接查看吗?
可以，用cat,less等命令均可查看
- Sam/Bam文件分为两部分，一部分以@开头的是什么序列，另一部分是什么序列？
以@开头的是头部区域，记录比对的整体信息；另一部分为主体区，记录详细的比对信息，至少11列
- Sam文件除头文件，以什么符号分割文本的，比对部分信息一行是多少列？你能用awk计算列数吗？
分隔符为tab，比对部分的列数必须有的是11列，剩余的列数不一定；用awk计算列数参考第13题答案
- Sam/Bam文件的@SQ开头的行是什么？你能生成一个文本，两列，一列是参考基因组染色体/sca id，一列是长度吗？
是参考基因组序列信息
- Sam文件的比对位置是从1还是0开始的？
从1开始
- 常见CIGAR 字符串各字母代表的意义
“M”表示 match或 mismatch；
“I”表示 insert
“D”表示 deletion
“N”表示 skipped（跳过这段区域）
“S”表示 soft clipping（被剪切的序列存在于序列中）
“H”表示 hard clipping（被剪切的序列不存在于序列中）
“P”表示 padding
“=”表示 match
“X”表示 mismatch（错配，位置是一一对应的）
- 比对质量的数字是哪一列？越大比对质量越好还是越小越好？
第五列，越大越好，60最好
- Sam文件的flag是第几列，数字代表什么意义？是怎么计算来的？
第二列，代表比对的情况，具体计算看表吧
- Sam文件怎么转bam文件？用什么指令？为什么要转换?
命令参考samtools view -bS test.sam > test.bam；Bam文件相对于Sam文件进一步压缩，同样的空间可以储存更多数据
- Bam文件查看用什么指令？如果需要查看头文件需要增加参数？
指令： samtools view; 查看头文件指令：samtools view test.bam | grep "^@"
- Bam文件为什么要排序？排序后的bam和未排序的bam头文件和chr position列有什么区别？
这个不是很确定···
- Bam文件建索引的指令是什么？
samtools index test.bam test.bam.bai
- Bam文件可视化用什么工具？查看时需要建立索引吗？
可视化可以用IGV吧
- Bam文件统计reads比对情况用samtools的哪一个子命令？
samtools tview
- SE测序和PE测序的所比对后得到的sam文件的区别在哪里？
根据对SAM文件格式的理解，Flag值会不一样；其他的不太确定···
- Bam能用gzip再压缩吗？
not sure about this
- Sam文件通常由哪些比对软件得到的？
Bowtie2啊，Hisat啊什么的
- Sam/Bam文件能转成fastq文件吗？
当然可以
- 有时候不能通过文件名的后缀来区别是sam还是bam，你是怎么区别的？
Sam文件可以直接用cat，less命令查看，Bam为二进制文件，需要用samtools view命令查看
