1. fasta和fastq

1.1. fasta：序列

1. 以 > 开头
2. gi|gi号|来源标识|序列标识（接收号/名称等）
3. 序列描述
4. 序列中允许空格，换行，空行，直到下一个 > 表示序列结束

1.2. fastq：序列+质量

1. 第一行以 @ 开头，唯一的序列ID标识符，可选的序列描述内容，用空格分开
2. 第二行：序列字符（核酸/氨基酸）
3. 第三行：“+”开头，后面空着或者加上第一行@之后的内容
4. 第四行：碱基质量字符，每个字符对应第二行相应位置碱基或氨基酸的质量，可以按一定规则转换为碱基质量得分，进而反映该碱基的错误率，这一行字符数与第二行字符数必须相等

1.3. 处理

• fastqc test.fq 即可得到html格式的报告，可以可视化地查看碱基是否平衡，测序质量如何等等

1.4. 作业

• 练习链接

1. echo $(wc -l reads_1.fq|awk '{print $1}')/4|bc
2. awk '/^@/{print $0}' reads_1.fq或awk '{if(NR%4==1)print}' read_1.fq
3. cat reads_1.fq |paste - - - - |cut -f2
4. cat reads_1.fq |paste - - - - |cut -f3
5. cat reads_1.fq |paste - - - - |cut -f4
6. cat reads_1.fq |paste - - - - |cut -f2|grep 'N'|wc
7. cat reads_1.fq |paste - - - - |cut -f2|wc -m
   1. 用上面这种方法会多出10000个回车符，所以要么手动减去10000，要么用下面的方法
   2. awk '{if(NR%4==2)print length}' read_1.fq|paste -sd +|bc
   3. cat reads_1.fq |paste - - - - |cut -f2|grep -o [ATCGN]|wc
8. cat reads_1.fq |paste - - - - |cut -f2|grep -o 'N'|wc
9. cat reads_1.fq |paste - - - - |cut -f4|grep -o '5'|wc
10. cat reads_1.fq |paste - - - - |cut -f4|grep -o '?'|wc
11. cat reads_1.fq |paste - - - - |cut -f2|grep '^[ATCGatcg]'|wc
    1. 理解错了题目意思，正解：awk '{if(NR%4==2)print length}' read_1.fq|cut -c1|sort|uniq -c
12. cat reads_1.fq |paste - - - - |cut -f1,2|tr '\t' '\n'|tr '@' '>'>reads_1.fa
13. echo $(wc -l reads_1.fa|awk '{print $1}')/2|bc
14. cat reads_1.fa |paste - - |cut -f2|grep -o [CG]|wc
15. cat reads_1.fa |tr -d 'N'
16. cat reads_1.fa |grep -v 'N'
17. cat reads_1.fa |paste - - |cut -f2|awk '{if(length>=65)print}'
    1. 这样会把第一行的序列信息去掉，用下面的方法更好
    2. cat reads_1.fa |paste - - |awk '{if(length($2)>=65)print}'
18. awk '{if(NR%2==0){print substr($0,6,length($0)-10)}}' reads_1.fa 这个做法的缺点同样是丢失了序列信息
19. cat reads_1.fa |paste - - |awk '{if(length($2)<=125)print}'
20. 直接用fastqc

测序质量.png

2. sam和bam

2.1. 标头注释部分

• 标头以@开头，用不同的tag表示不同的信息，主要有：
  • @HQ，说明符合标准的版本、对比序列的排列顺序
  • @SQ，参考序列说明
  • @PG，使用的比对程序说明
  • LN，参考序列的长度

2.2. 比对结果部分

1. 每一行表示一个read的比对信息
2. 每行包括11个必须字段和1个可选字段，字段之间用tag分割

2.3. 作业

• 练习链接

1. cat tmp.sam |grep -v '@'|cut -f1|uniq|wc
2. cat tmp.sam |grep -v '@'|cut -f2|sort|uniq -c|sort -k1,1 -n
3. cat tmp.sam |grep -v '^@' |awk '{if($6=="*")print$10}'|head -3

cat tmp.sam |grep -v '^@' |awk '{if($6=="*")print$1}'|sort|uniq -c|grep -w 1
cat tmp.sam |grep -v '^@' |awk '{if($6=="*")print$1}'|sort|uniq -c|grep -w 2

5. cat tmp.sam |grep -v '^@' |awk '{if($5>30)print}'
6. cat tmp.sam |grep -v '^@' |awk '{if($6!="*")print$6}'|grep '[IDNSHPX]'
7. cat tmp.sam |grep -v '^@' |awk '{if($9>1250||$9<-1250)print}'
8. cat tmp.sam |grep -v '^@' |cut -f3 |sort|uniq -c
9. cat tmp.sam |grep -v '^@' |awk '{if($10 ~/N/) print $0}'
10. awk '{if($10 ~/N/) print}' tmp.sam |awk '{if($6 !~ /[IDNSHPX*]/) print}'
11. grep '^@' tmp.sam |wc
12. grep -v '^@' tmp.sam |cut -f 12-1000|head -6
13. grep -v '^@' tmp.sam |cut -f 12-1000|awk '{print NF}' |sort |uniq -c
14. grep '^@SQ' tmp.sam |awk '{print $3}'|cut -d: -f2
15. awk '{if($6 ~/I/)print}' tmp.sam
16. awk '{if($6 ~/D/)print}' tmp.sam
17. awk '{if($4>5013 && $4<50130) print}' tmp.sam
18. cat tmp.sam |grep -v '^@' |awk '{if($4!=0)print}'|sort -nk4,4|less -SN
19. grep "102M3D11M" tmp.sam |awk '{print length($10)}'

2.4. 其它概念题

1. 高通量测序数据分析中，序列与参考序列的比对后得到的标准文件为什么文件？

A：sam文件

2. sam文件是一种比对后的文件，能直接查看吗？

A：可以

3. Sam/Bam文件分为两部分，一部分以@开头的是什么序列，另一部分是什么序列？

A：分别是标头注释部分和比对结果部分

4. Sam文件除头文件，以什么符号分割文本的，比对部分信息一行是多少列？你能用awk计算列数吗？

A：以tab分割，必须字段11列，可选字段不一定，可以用awk '{print NF}'

5. Sam/Bam文件的@SQ开头的行是什么？你能生成一个文本，两列，一列是参考基因组染色体/sca id，一列是长度吗？

A：参考染色体的信息，cat tmp.sam|grep -w ^@SQ|cut -f2-3 > chrom.txt

6. Sam文件的比对位置是从1还是0开始的？

A：从1开始

7. 常见CIGAR字符串各字母代表的意义

A：

• M：alignment match (can be a sequence match or mismatch)表示read可mapping到第三列的序列上，则read的碱基序列与第三列的序列碱基相同，表示正常的mapping结果，M表示完全匹配，但是无论reads与序列的正确匹配或是错误匹配该位置都显示为MI：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”表示错误匹配到序列上

8. 比对质量的数字是哪一列？越大比对质量越好还是越小越好？

A：第五列，越大越好。

9. Sam文件的flag是第几列，数字代表什么意义？是怎么计算来的？

A：第二列

• 1（1）该read是成对的paired reads中的一个
• 2（10）paired reads中每个都正确比对到参考序列上
• 4（100）该read没比对到参考序列上
• 8（1000）与该read成对的matepair read没有比对到参考序列上
• 16（10000）该read其反向互补序列能够比对到参考序列
• 32（100000）与该read成对的matepair read其反向互补序列能够比对到参考序列
• 64（1000000）在paired reads中，该read是与参考序列比对的第一条
• 128（10000000）在paired reads中，该read是与参考序列比对的第二条
• 256（100000000）该read是次优的比对结果
• 512（1000000000）该read没有通过质量控制
• 1024（10000000000）由于PCR或测序错误产生的重复reads
• 2048（100000000000）补充匹配的read

多个数字相加得到

10. Sam文件怎么转bam文件？用什么指令？为什么要转换？

A：samtools view -bS tmp.sam >tmp.bam 因为bam占用的磁盘空间比sam文本文件小；利用bam二进制文件的运算速度快

11. Bam文件查看用什么指令？如果需要查看头文件需要增加参数？

A：samtools view tmp.bam，加上选项-H

12. Bam文件为什么要排序？排序后的bam和未排序的bam头文件和chr position列有什么区别？

A：排序后所有reads按照其在参考染色体上的位置排好

13. Bam文件建索引的指令是什么？

A：samtools index <in.bam> [out.index]

14. Bam文件可视化用什么工具？查看时需要建立索引吗？

A：igv 需要

15. Bam文件统计reads比对情况用samtools的哪一个子命令？

A：samtools idxstats tmp.bam

16. SE测序和PE测序的所比对后得到的sam文件的区别在哪里？

A：qname第一列，双端有重复现象，单端没有重复现象

17. Bam能用gzip再压缩吗？

A：不能，虽然可以执行tar -zcvf tmp.tar.gz tmp.bam命令，但随后用ls -lh查看可以发现文件大小没有变化

结果示意.png

18. Sam文件通常由哪些比对软件得到的？

A：bowtie2，bwa

19. Sam/Bam文件能转成fastq文件吗？

A：能

samtools view -bS aln.sam > aln.bam
samtools view -h -o aln.sam aln.bam
bam2fastq --aligned input.bam -o output.fq

20. 有时候不能通过文件名的后缀来区别是sam还是bam，你是怎么区别的？

A：能否直接cat查看

3. gff和gtf

• gff主要用来注释基因组
• gtf主要用来注释基因

3.1. gtf

• 每一列的含义
  1. seqid 参考序列ID
  2. source 注释的来源 未知则用 . 代替。一般指明产生此文件的软件或方法
  3. type 类型 此处名词相对自由，建议使用符合SO管理的名称，如gene，repeat_region，exon，CDS等
  4. start 开始位点 从1开始计数
  5. end 结束位点 对于一些可以量化的属性，可以在此设置一个数值以表示程度的不同。如果为空，用 . 代替。
  6. score 得分
  7. strand 正/负链
  8. phase 步进 对于编码蛋白质的CDS，本列指定下一个密码子开始的位置，可以是0、1或2，表示到达下一个密码子要跳过的碱基个数
  9. attributes 属性 格式为标签=值（tag=value）不同属性之间以分号相隔

3.2. gff

• 每一列的含义
  1. seqname 序列的名字 通常格式染色体ID或是contig ID
  2. source 注释的来源 通常是预测软件名或是公共数据库
  3. start 开始位点 从1开始计数
  4. end 结束位点
  5. feature 基因结构 CDS，start_condon，stop_condon是一定要含有的类型
  6. score 得分 对该类型存在性和其坐标的可信度，非必须，可以用点 . 代替
  7. strand 正/负链 链的正向与负向，分别用+/-表示
  8. frame 密码子偏移 可以是0，1或2
  9. attributes 属性
     1. gene_id value：表示转录本在基因组上的基因座的唯一的ID。
     2. gene_id与value值用空格分开，如果值为空，则表示没有对应的基因。
     3. transcript_id value：预测的转录本的唯一ID
     4. transcript_id与value值用空格分开，空表示没有转录本。

3.3. 关系

• gtf2和gff3的内容很相似，区别只在其中的两列

• 两种文件格式之间的转换：使用Cufflinks里的工具“gffread”

#gff2gtf
gffread my.gff3 -T -o my.gtf
#gtf2gff
gffread merged.gtf -o- >merged.gff3

4. vcf

• VCF格式，Variant Call Format 变异判读文件格式
• 分为两部分内容：以“#”开头的注释部分；没有“#”开头的主体部分
• VCF文件主题部分的结构
  1. CHROM：表示变异位点是在哪个contig里call出来的，如果是人类全基因组的话那就是chr1…chr22, chrX, Y, M
  2. POS： 变异位点相对于参考基因组所在的位置，如果是indel，就是第一个碱基所在的位置
  3. ID： 如果call出来的SNP存在于dbsnp数据库里，就会显示相应的dbsnp里的rs编号, 若没有，则用 . 表示其为一个novel variant。
  4. REF： 在这个变异位点处，参考基因组中所对应的碱基
  5. ALT： 在这个变异位点处，研究对象基因组中所对应的碱基
  6. QUAL：Phred格式(Phred_scaled)的质量值，表示在该位点存在variant的可能性；该值越高，则variant的可能性越大；计算方法：Phred值 = -10 * l/g (1-p) p为variant存在的概率; 通过计算公式可以看出值为10的表示错误概率为0.1，该位点为variant的概率为90%。
  7. FILTER：使用上一个QUAL值来进行过滤的话，是不够的。GATK能使用其它的方法来进行过滤，过滤结果中通过则该值为'PASS'；若variant不可靠，则该项不为'PASS'或'.'。
  8. INFO
  9. FORMAT 和 testxxx：这两行合起来提供了'testxxx'这个sample的基因型的信息。 testxxx 代表这该名称的样品，是由BAM文件中的@RG下的 SM 标签决定的。
     • GT：样品的基因型（genotype）。两个数字中间用'/'分开，这两个数字表示双倍体的sample的基因型。0 表示样品中有ref的allele；1 表示样品中variant的allele； 2表示有第二个variant的allele。因此：0/0 表示sample中该位点为纯合的，和ref一致；0/1 表示sample中该位点为杂合的，有ref和variant两个基因型；1/1 表示sample中该位点为纯合的，和variant一致。
     • AD：对应两个以逗号隔开的值，这两个值分别表示覆盖到REF和ALT碱基的reads数，相当于支持REF和支持ALT的测序深度。
     • DP：覆盖到这个位点的总的reads数量，相当于这个位点的深度（并不是多有的reads数量，而是大概一定质量值要求的reads数）。
     • PL：指定的三种基因型的质量值(provieds the likelihoods of the given genotypes)。这三种指定的基因型为(0/0,0/1,1/1)，这三种基因型的概率总和为1。和之前不一致，该值越大，表明为该种基因型的可能性越小。 Phred值 = -10 * lg § p为基因型存在的概率=10^（-Phred值/10）。
       INFO
     • AC：表示该Allele的数目，Allele数目为1表示双倍体的样本在该位点只有1个等位基因发生了突变
     • AF：表示Allele的频率，Allele频率为0.5表示双倍体的样本在该位点只有50%的等位基因发生了突变
     • AN：表示Allele的总数目,即：对于1个diploid sample而言：则基因型 0/1 表示sample为杂合子，Allele数为1(双倍体的sample在该位点只有1个等位基因发生了突变)，Allele的频率为0.5(双倍体的 sample在该位点只有50%的等位基因发生了突变)，总的Allele为2； 基因型 1/1 则表示sample为纯合的，Allele数为2，Allele的频率为1，总的Allele为2。
     • DP：样本在这个位置的reads覆盖度，是一些reads被过滤掉后的覆盖度（跟上面提到的DP类似）
     • FS：使用Fisher's精确检验来检测strand bias而得到的Fhred格式的p值，值越小越好
     • MQ：表示覆盖序列质量的均方值RMS Mapping Quality

4.1. 作业

• 练习链接

1. cat tmp.vcf|grep -v "^#"|grep INDEL``cat tmp.vcf|grep -v "^#"|grep -v INDEL

(test) [root@instance-ptelbf6dreads]# cat tmp.vcf|grep -v "^#"|grep INDEL|cut -f8|cut -d';' -f4|cut -d'=' -f2|paste -sd +|bc
2831
(test) [root@instance-ptelbf6dreads]# cat tmp.vcf|grep -v "^#"|grep INDEL|wc -l
70
(test) [root@instance-ptelbf6dreads]# echo 2831/70|bc
40
#以上方法可能略笨，而且只能得到整数，下面附上别人写的更好的方法
(test) [root@instance-ptelbf6dreads]# cat tmp.vcf|grep -v "^#"|grep INDEL|cut -f8|cut -d";" -f4|cut -d"=" -f2|awk '{ sum += $0; } END { print sum/NR }'
40.4429
(test) [root@instance-ptelbf6dreads]# cat tmp.vcf|grep -v "^#"|grep -v INDEL|cut -f8|cut -d";" -f1|cut -d"=" -f2|awk '{ sum += $0; } END { print sum/NR }'
37

cat tmp.vcf|grep -v "^#"|grep INDEL|awk '{if(length($4)<length($5))print}' # insertion
cat tmp.vcf|grep -v "^#"|grep INDEL|awk '{if(length($4)>length($5))print}' # deletion

4. cat tmp.vcf|grep -v "^#"|grep -v INDEL|awk '{print $4"-->"$5}'|sort|uniq -c|sort -k1,1nr SNP行的第4、第5行分别表示参考基因和研究对象的碱基类型
5. cat tmp.vcf|grep -v "^#"|cut -f10|cut -d: -f1|sort|uniq -c 缺陷是没有把整条取出来
6. ``暂略
7. cat tmp.vcf|grep -v "^#"|awk '{if($6>30)print}'
8. ``暂略
9. ``暂略
10. IGV暂略

4.2. 其他思考题

1. vcf的全称是什么？是用来记录什么信息的标准格式的文本？

A：VCF文件全称为Variant Call Format，表示基因组的变异信息。

2. 一般选用哪个指令查看vcf文件，为什么不用vim?

A：一般选用less -SN，因为vim没有单行显示，看得不清楚。

3. vcf文件以'##'开头的是什么信息？请认真查看这些信息。'#'开头的是什么信息？

A：以##开头的是mate-information line，以##开头，格式为key=value。 fileformat是必须的字段，表明VCF格式的版本，其他行主要用来描述INFO, FORMAT, FILTER等字段的具体含义；以#开头的是heaher line，只有一行，行内以\t分隔，是下方具体信息的表头。

4. vcf文件除头信息，每行有多少列，请详细叙述每行的含义！请准确记忆。

A：至少9列。1）CHROM: 染色体名字；2）POS: 染色体的位置，起始位置为1；3）ID: 变异位点在数据库中的ID，如果是dbsnp数据库，推荐使用rs号，如果没有ID，用点号表示缺失值；4）REF: 参考基因组上的碱基；5）ALT: 变异之后的碱基；6）QUAL: 变异位点的质量，质量值越高，为真实的变异位点的概率越大；7）FILTER: 过滤信息，PASS代表通过了过滤；对于过滤失败的位点，会给出对应的过滤失败的原因，具体的含义可以查看mate-information line 中对FILTER字段的描述；8）INFO：额外的信息，具体的含义可以查看mate-information line 中对INFO 字段的描述，看上去是一列，但其中的内容可以无限扩增；9）Format：表头的##FORMAT就是对第九列的解释，主要包括某一个特定位点基因型、测序深度的表述。

5. 理解format列和样本列的对应关系以及GT AD DP的含义。

A：GT：genotype；AD：Allele Depth，样本中每一种allel的reads覆盖度，在二倍体中是用逗号分隔的两个数，前面对应ref，后面对应variant；DP：Raw read depth

6. vcf文件第三列如果不是'.'，出现的rs号的id是什么？

A：变异位点名称，对应dbSNP数据库中的ID

7.** vcf文件的ref，alt列和样本列的0/1 1/1 或者1/2的联系？**

A：0代表样本中ref的allel；1代表样本variant的allel；2表示有第二个variant的allel；0/0表示样本中该位点纯合，与ref一致；0/1表示样本中该位点杂合，有ref和variant两个基因型；1/1表示样本中位点纯合，与variant一致

8. vcf文件一般用什么软件生成？请至少说出两个软件。请注意不同软件生成的vcf格式的稍有不同的地方。

A：GATK和Samtools

9. Vcf文件一般都比较大，压缩后的.gz文件用什么指令直接查看而不用解压后查看？

A：zless zcat

10. 了解gvcf是什么格式，gvcf全称是什么？他与vcf有什么前后联系？

A：GVCF代表Genomic VCF，GVCF是一种VCF，因此基本格式规范与常规VCF相同，但基因组VCF包含额外信息。术语GVCF有时仅用于描述包含基因组中每个位置（或感兴趣区间）的记录的VCF，无论是否在该位点检测到变体（例如由UnifiedGenotyper产生的VCF --output_mode EMIT_ALL_SITES）。HaplotypeCaller 3.x生成的GVCF包含以非常特定的方式格式化的其他信息。

11. 把alt列出现','的行提取出来

A：cat tmp.vcf |grep -v "^#"|awk '{if($5 ~",")print}'

12. 请将chrid、postion、ref、alt、format、样本列切割出来生成一个文本

A：cat tmp.vcf |grep -v "^#"|awk '{print $1,$2,$4,$5,$9}'> tmp.new.txt

13. 将一个含snp，indel信息的vcf拆成一个只含snp，一个只含indel信息的2个vcf文件。可借鉴软件

A：

cat tmp.vcf |grep -v "^#"|awk '{if($8~"INDEL")print}'>indel.vcf
cat tmp.vcf |grep -v "^#"|awk '{if($8!~"INDEL")print}'>snp.vcf

14. 用指令操作indel的vcf文件，提取indel长度>4的变异行数，存成一个文本。

A：cat indel.vcf |cut -f8|cut -d';' -f2无法提取完整信息

15. 用vcftools过滤vcf文件，如maf设置成0.05， depth设置成5-20，统计过滤前后的变异位点的总个数

16. 利用vcftools提取每个样本每一个位点的变异信息和深度信息，生成一个矩阵的文件，至少含chrid，postion，sample_DP，sample_GT四项信息

vcftools暂略

17. 提取出变异位点上样本有纯和突变的行数

A：grep 1/1 snp.vcf |wc

18. 统计一下1号染色体上的变异总个数。

A：示例文件的序列全部在一条染色体上。

19. 提取一下BRCA1基因上发生变异的行数，如果是人类wes变异结果文件

暂略

20. 统计vcf文件各样本的缺失率，如果是多个样本的群体call结果。

暂略

5. 小技巧

• find . -name *.sam 搜索sam格式文件

6. 本文引用源

1. https://blog.csdn.net/sunchengquan/article/details/86677629
2. https://www.jianshu.com/p/396280f55558

最后，向大家隆重推荐生信技能树的一系列干货！

1. 生信技能树全球公益巡讲：https://mp.weixin.qq.com/s/E9ykuIbc-2Ja9HOY0bn_6g
2. B站公益74小时生信工程师教学视频合辑https://mp.weixin.qq.com/s/IyFK7l_WBAiUgqQi8O7Hxw
3. 招学徒https://mp.weixin.qq.com/s/KgbilzXnFjbKKunuw7NVfw