刘小泽写于19.6.3
继续整理直播基因组优质内容，这一次介绍基因组分析两种基本思路
基因组的那些事儿开篇：https://www.jianshu.com/p/bf871522ea20
思路一：科研分析
测完了基因组，会得到大量的数据，最有价值的分析就是通过与参考基因组的比对，来找出不一样的地方，得到VCF文件（储存序列变异的文本文件格式，tab分割）。另外，并非所有与参考基因组不一样的地方都是有意义的，还需要根据一系列的数据库来注释找到variation 。
通常一个人的全基因组测序数据可以挖掘到四百万个SNVs(和参考基因组不一样的单碱基位点)，还有五十万的indels(insertions or deletions），这些信息都被存储在VCF文件中。一开始我们只能知道20号染色体上14370这个位点的G变成了A，仅此而已。但VCF其实很多，比如，我们发现G变成了A，那么这个位点是在某个基因上面吗？如果是，在基因的外显子还是内含子？它的变异有没有改变这个基因的功能呢？有没有影响它的转录和翻译呢？还有没有其他正常人也是这个位点变异呢？如果有，是哪些人种呢？有没有癌症病人也发现了这个变异呢？如果有，是什么癌症呢？～是吧，蕴藏的信息量巨大，而这些信息的挖掘需要依靠变异位点注释数据库
这些变异位点注释数据库主要包括：
dbsnp（ncbi提供的最权威）:https://www.ncbi.nlm.nih.gov/projects/SNP/
ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh38p7/
ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh37p13/
ExAC（Exome Aggregation Consortium）[broadinstitute提供的外显子联盟]:http://exac.broadinstitute.org/
ftp://ftp.broadinstitute.org/pub/ExAC_release/current
Cosmic癌症突变信息集【需要学术邮箱注册】
TCGA:最大的癌症基因信息的数据库,搜集的是TCGA计划里面总结的各个癌症somatic mutation【样本的变异文件和它有交集，就不是什么好事】
1000genomes【千人基因组计划】 ：包括5个大人种，共25个小人种的基因型数据，结果会得到一定程度上的祖缘分析结果
ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/
ESP（Exome Variant Server）
http://evs.gs.washington.edu/EVS/国家级数据库 ：
UK10K 英国一万人基因组：https://www.uk10k.org/
GoNL(Genome of the Netherlands) 荷兰人：http://www.nlgenome.nl/
https://molgenis26.target.rug.nl/downloads/gonl_public/variants/release6.1/
SSMP（Singapore Sequencing Malay Project） 马来西亚
SSIP（Singapore Sequencing Indian Project） 印度
数据下载、质控、过滤QC-》比对alignment（bwa）-〉sam转bam并排序sort-》标记PCR重复，并去除-〉产生需要重排的坐标记录-》根据重排记录文件把比对结果重新比对-〉最终的bam文件转为mpileup文件-》variation calling（bcftools+freebayes+gatk+varscan）-》annotation-statistics/visualization-》snv,indel,cnv,sv
其中将测序数据与参考基因组进行比对（alignment）环节对计算资源要求较高，比对的速度取决于cpu的配置以及软件，一般8线程一秒钟也才能比对1000条双端测序reads。
质控用的软件：fastqc、multiqc【一次同时生成多个数据质量报告】、fastp
过滤用的软件：
trim galore【和fastqc出自一家，可以和fastqc结合使用，用来清洗原始数据】
trimmomatic【专门清洗illumina测序数据】
fastp：国内出的软件，简单易懂【仅仅扫描 FASTQ 文件一次，就完成比FASTQC + cutadapt + Trimmomatic 这三个软件加起来还多很多的功能，而且速度上比仅仅使用 Trimmomatic 一个软件还要快 3 倍左右；支持多线程】
为什么要做质控、过滤这个事？
常用的二代测序基本都是illumina的边合成边测序技术。怎么合成？靠的是DNA聚合酶是碱基链延伸，我们大家长时间工作都会疲惫，这个酶也一样，合成超过一定长度的链错误率就会升高。不管怎样，结果出来了，你拿到数据怎么知道测的合不合格呢？
这就需要质控。另外还需要了解，不同测序平台都有测序错误偏向性（比如illumina平台对reads的5‘端的前几个随机引物碱基有一定的偏好性，因此碱基分布图中前端波动较大）
有问题就解决～一般需要去除低质量碱基和接头
比对、排序 ：
BWA（Burrow-Wheeler Aligner）：李恒开发的使用最广的NGS数据比对软件，目前更新到0.7.17
为什么要做比对这个事？
本来好端端的短序列read都十分有序的躺着基因组中，有一天你拿来测基因组，但经过DNA建库和测序后，原有的平静被打破，不同read之间的顺序关系也混乱不堪，因此你没办法判断两条read谁前谁后，这时你只能把他们当成初始read1、初始read2；
但是比对软件给了你希望，他可以将一大批read一个一个与参考基因组比较（这里的比对策略一般有overlap、k-mer），再按顺序排列，整齐划一的队伍才好管理！
说到比对，不知你是否想到了我们最常用的ncbi blast？你每次提交一个序列，结果就会给你一堆对应的其他序列，看他们相似性。这就是比对的过程：你提交的是一个字符串，比对库就是大的字符串合集，因此，比对就是搜索最相似的字符串的过程。但是利用blast的动态规划算法，比对基因组的reads，那是一个十分缓慢的过程（你可以回想一下提交一条序列去blast需要多久，而一般基因组测序会有6亿条150bp左右短序列），因此，这里的比对也要特殊优化。
BWA为什么是一个优秀的比对软件呢？因为他采用了BW（Burrow-Wheeler）这一套压缩算法，先对参考基因组进行构建索引。这里可以假设一下，你拿着10个书名，去图书馆借书，到了之后，有两种方法可以拿到这10本书：一是从第一排开始一本一本找；二是根据图书馆图书查找系统，准确定位。我们这里的参考基因组就是一个图书馆，我们要先为他打造一个查找系统，就是index索引。BWA就是使用（BW）块排序压缩的算法完成这个工作
Samtools：同样由李恒开发，用作数据格式转换、比对、排序
BWA得到的bam没有顺序么？为什么要排序？
fastq文件中reads在基因组上随机分布，也就是说fastq文件中的reads之间是没有顺序的，而上一步比对仅仅是是将不同fastq文件中的reads定位到参考基因组，定位完就输出bam文件了，并没有自动识别reads的前后位置并重新排序。因此，在初步得到的bam文件中，每一条比对记录都是没有先后顺序的。
但是呢，后来不管是要IGV可视化还是要去重复，都需要bam中的记录有顺序，还是那句话“整齐划一的才好管理！”
寻找变异位点：bcftools、GATK、freebayes、VarScan
注释用到的软件：VEP，ANNOVAR, SnpEFF
思路二：临床流程
我们得到了许多变异位点，那么这些不同的地方是什么意思呢？这就需要和公共疾病数据库比对进行整理，注释那些不一样的地方。OMIN，clinVAR，HGMD，GWAS。大部分疾病评估是依据GWAS数据库对变异位点进行注释从而评估个体化疾病风险；用药建议是根据PharmaGKB网站；遗传病风险则是HGMD数据库。
变异与突变是不同的！
变异只是产生人种多样性的原因， 有了变异我们才能不一样；而那些跟疾病相关的变异，通常才叫做突变。一般，在正常人的数据库里面出现了5%的变异就可以认为没什么大的危害
变异分成4种，即snv、indel、cnv、sv，大部分情况下只能分析到SNV。【另外3个要么不准确，要么有点难度】
欢迎关注我们的公众号～_～
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到jieandze1314@gmail.com
