【好久没记录了~最近在准备复试的时候发现了远古时期自己做的大作业，虽然漏洞百出（那时候居然会把个体叫做1位个体），但maybe是我大学期间最用心的大作业了。谁让邵老师魅力无穷呢】

遗传学是破译不同人群表型多样性的关键。近年来，基于microarray和第二代测序技术 ，Hapmap 和千人基因组计划取得了重大的成就。韩国人口是 HapMap 和千人基因组计划中未包括的亚洲人口， 目前还没有一个综合的韩国群体的遗传结构来解释韩国人群特有的遗传特征。本文旨在利用全基因组测序数据来描述韩国人群的遗传特征。 通过与 HapMap 和 千人基因组计划 检测到的遗传变异进行比较分析，我们确定了只存在于韩国人群中的 snv ，并探讨了它们与功能、通路 、疾病和药物方面的关系 。这些发现加深了我们对韩国人群遗传学和进化的理解，并有望促进韩国人群的 个性化医疗 。

1、 CAMDA 获取数据（ http://dokuwiki.bioinf.jku.at/doku.php/start ）。数据包括全基因组测序原始读段，BWA比对结果，用 SAMTools 做的 SNV calling 结果。从中选取 3 5个韩国血统人群样本用作研究。 2、 由全基因组测序数据得到的 9 个韩国人的 SNVs hg19格式（http://tiara.gmi.ac.kr/download） 3、 hg19 参考序列（ftp://hgdownload.cse.ucsc.edu /goldenPath/hg19/bigZips/chrom FaMasked.tar.gz）

本文使用两种SNV calling的方法分析35位韩国人的全基因组测序的原始读段。 第一种用BWA（0.5.9）与人类基因组（hg19）进行读段比对，然后用SAMtools进行SNV calling；第二种用SOAP2（2.21）（bowtie2）与人类基因组（hg19）进行读段比对并且用SOAPsnp（GATK）做SNV calling。将用两种方法的得到的 SNV 根据基因位置合并在一起，然后进行比较。提取 那些 用两种方法均能检测到的被认为是高质量的韩国人群 SNV 。 接下来，比较韩国 SNV 和在 1KGP 、 HapMap 中 能 检测到的其他人群的 SNV并将韩国人群 SNV 分为两类：仅在韩国人群的 SNV Korean与韩国人群和其他人群共享的 SNV 。最后 根据注释，确定仅在韩国 人群 的 非同义 SNV ，然后鉴定涉及的基因以进行 GO和 KEGG 富集分析，并探讨它们与疾病和药物的关联 。

由于全基因组数量过于庞大，于是只使用了第一位韩国人的全基因组进行本 次研究。即 KPGP_ 00001。

根据作者提供的网址，去 CAMDA 下载，发现已经失效了！于是又去 KPGP 千人基因组计划官网查询，发现全是韩文！接着又去http://opengenome.net/index.php/Main_Page ，发现自己被 forbbidden 了。最后非常绝望，只能去万能的百度上搜索，在简书中遇见了KPGP 基因组！！

注： wget 的速度非常非常非常非常 慢！但因为这次下载的数据比较特殊，下 次如果下载 UCSC 、GEO 上数据，可以使用 sratoolkit ，或者强大的迅雷下载

"Two sets of alignment results were generated for the raw reads of the 35 Korean samples” 本文未对测序读段进行质控，故我也不用质控啦。

一、 使用BWA做比对 “The first set was provided by KPGP and was generated by using BWA (version 0.5.9)to map raw reads to the human genome (hg19) with 45bp seed sequence allowed (see Additional file 1 for details).”

二、使用bowtie2做比对（原文使用的是soap2，但因华大已经两年未更新，且相应下游软件SOAPsnp已经不维护了，改用bowtie2作比较）

一、bwa结果进行SNP calling

删除pcr前后比对结果如下

分别是96.44%与94.04%

查看vcf文件

注：vcf格式解释 CHROM和POS：变异位点所在的染色体名称和位置，从１开始计数，如果是INDEL的话，位置是该INDEL第一个碱基的位置。 ID：variant的id。比如call出来的SNP在dbSNP数据库中存在，这里就会显示相应的rs号（当然前提是已经和dbSNP数据库做了比较）。 REF和ALT：参考序列的碱基和突变后的碱基。如果有多种不同于参考序列的基因型，在ALT列使用“，”隔开。如变异位点在参考基因组上的碱基为“G”，样品上突变后的基因型为“A”，则REF列为“G”，ALT列为“A”；如果突变后的碱基有多个如A和C，则ALT可以表示为“A,C”。这里需要注意ALT是针对这个变异位点而言，不针对特定样品。 QUAL：Phred格式(Phred_scaled)的质量值，表示在该位点存在variant的可能性，值越高，则variant的可能性越大。计算方法：Phred值= -10 * log (1-p)， p为variant存在的概率。通过计算公式可以看出值为10的表示错误概率为0.1，该位点为variant的概率为90%。 FILTER：理想情况下，QUAL这个值应该是用所有的错误模型算出来的，这个值就可以代表正确的变异位点了，但是事实是做不到的。因此，还需要对原始变异位点做进一步的过滤。无论你用什么方法对变异位点进行过滤，过滤完了之后，在FILTER一栏都会留下过滤记录，如果是通过了过滤标准，那么这些通过标准的好的变异位点的FILTER一栏就会注释一个PASS，如果没有通过过滤，就会在FILTER这一栏提示除了PASS的其他信息。如果这一栏是一个“.”的话，就说明没有进行过任何过滤。 INFO：这一列是variant的详细信息，格式以tag=value形式记录，而tag的说明一般包含在文件开头的注释说明部分。 FORMAT和NA00001（NA00002…）：FORMAT这列规定了后边样品每列的格式，NA00001（NA00002…）等各列是对应每个样品在这个variant的信息。我们如果要看每个样品的基因型信息，就需要看这几列了。

计算SNP总数

得出3,208,582个SNP

二、 bowtie2结果进行SNP calling（使用GATK） 附：GATK流程 注：GATK运行之前有一套复杂且必须的步骤要执行，不能偷懒

查找出3,413,019个SNP

BWA准确率似乎高一些，然后bowtie2速度快许多 GATK能识别到更多SNP位点，但是非常复杂。 样本数太少，不能得出有效结论

总结：不在千人基因组计划中的国外数据下载好难！GATK好麻烦，但应该比较准确！要记住bwa的index和参考序列要用统一的！！！