遗传学是破译不同人群表型多样性的关键。近年来，基于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 基因组！！

#下载了KPGP_志愿者的全基因组0001存在于data文件夹中 nohup wget ftp://biodisk.org/Release/KPGP/KPGP_Data_2013_Release_Candidate/WGS/K PGP 00001/KPGP 00001_L1_R1.fq.gz 1>/dev/null 2>&1 nohup wget ftp://biodisk.org/Release/KPGP/KPGP_Data_2013_Release_Candidate/WGS/K PGP 00001/KPGP 00001_L1_R2.fq.gz 1>/dev/null 2>&1 

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

"Two sets of alignment results were generated for the raw reads of the 35 Korean samples” 本文没有质量控制测序读段，所以我不需要质量控制。

#如果要用的话 fastp -i KPGP-00001_L1_R1.fq.gz -I out. KPGP-00001_L1_R1.fq.gz -I KPGP-00001_L1_R2.fq.gz -O out.KPGP-00001_L1_R2.fq.gz 

一、 使用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 llowed (see Additional file 1 for details).”

#之前没有建立过hg19索引，所以首先去UCSC下载hg19
wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.Masked.tar.gz
#注：没必要，老师目录里有建好的hg19索引
cd bwa
tar -zxf chromFa.tar.gz #解压参考基因组，显示所有染色体单个文件
cat *.fa >ref.fa #将解压后的.fa文件进行合并
bwa-build ref.fa #对hg19文件进行索引

#用bwa进行比对
bytlib load bwa.kit_0.7.15
bwa index ref.fa
bwa mem ref.fa ../data/KPGP-00001_L1_R1.fq.gz ../data/KPGP-00001_L1_R2.fq.gz > KPGP_00001.bwa.pe.sam
#注：BWA，bowtie2运用了BWT算法，所建index一定要和后面比对的基因组保持一致，不然就相当于白建了index（第一次就白建了）

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

#bowtie2建立索引
bowtie2 index /public/workspace/shaojf/Course/NGS/Reference/bowtie2_Index/GRCh38.bowtie2
#用bowtie2进行比对
nohup bowtie2 --phred -p 6 -x /public/workspace/shaojf/Course/NGS/Reference/bowtie2_Index/GRCh38.bowtie2 -1 KPGP_00001_L1_R1.fq.gz -2 KPGP_00001_L1_R2.fq.gz -S KPGP_00001.bowtie2.pe.sam

一、bwa结果进行SNP calling

#用view将bam转成sam
samtools view -b KPGP_00001.bwa.pe.bam -@ 2 KPGP_00001.bwa.pe.sam
#对bam文件排序
samtools sort -o KPGP_00001.bwa.pe.srt.bam -@ 2 KPGP_00001.bwa.pe.srt.rmdup.bam
#删除PCR重复
samtools rmdup -S KPGP_00001.bwa.pe.srt.bam KPGP_00001.bwa.pe.srt.rmdup.bam
#给原始转成bam的比对文件和删除pcr重复的bam文件建立索引，生成的索引文件后缀是bai
samtools index KPGP_00001.bwa.pe.srt.bam
samtools index KPGP_00001.bwa.pe.srt.rmdup.bam
#查看读段比对情况
nohup samtools flagstat KPGP_00001.bwa.pe.bam >rmdup.befor.flagstat 1>/dev/null 2>&1
nohup samtools flagstat KPGP_00001.bwa.pe.srt.rmdup.bam > rmdup.after.flagstat 1>/dev/null 2>&1

删除pcr前后比对结果如下 分别是96.44%与94.04%

#用samtools mpileup和bcftool call查看SNP位点（方法不太好，可以用bcftools mpileup直接替代）
samtools mpileup -f re.fa -o lab5.mpileup.vcf -uv lab4.bwa.pe.srt.rmdup.bam
bcftools call -mv --threads 10 KPGP_00001.mpileup.vcf -o KPGP_00001.raw.vcf 1>/dev/null 2>&1

查看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运行之前有一套复杂且必须的步骤要执行，不能偷懒

#将比对完的sam文件转成bam文件并排序
samtools view -b -o KPGP_00001.bowtie2.pe.bam -@ 10 KPGP_00001.bowtie2.pe.sam
nohup samtools sort -o KPGP_00001.bowtie2.pe.srt.bam -@ 10 KPGP_00001.bowtie2.pe.bam

#把sam文件转成bam
samtools view -b -o KPGP_00001.bwa.pe.bam -@ 2 KPGP_00001.bwa.pe.sam #编辑#ReadGroups信息
picard AddOrReplaceReadGroups I=KPGP_00001.bowtie2.pe.bam O= KPGP_00001.bowtie2.pe.rg2.bam ID=HCT116 LB=WXS PL=Illumina PU=Run SM=whole RGSM=20 #标记PCR重复之前先coordinate sort
picard SortSam I=KPGP_00001.bowtie2.pe.rg2.bam O=KPGP_00001.bowtie2.pe.rg2.srt.bam SORT_ORDER=coordinate #标记PCR重复
picard MarkDuplicates I=KPGP_00001.bowtie2.pe.rg2.srt.bam O=KPGP_00001.bowtie2.pe.rg2.md.bam M=KPGP_00001.bowtie2.marked_dup_metrics.txt ##执行BQSR
Time gatk BaseRecalibrator -R /public/workspace/shaojf/Course/NGS/Reference/bowtie2_Index/GRCh38.bowtie2 -I KPGP_00001.bowtie2.marked_dup_metrics.txt -O KPGP_00001.sorted.markdup.BQSR.bam #调用HaplotypeCaller输出样本VCF
Time gatk HaplotypeCaller -R /public/workspace/shaojf/Course/NGS/Reference/bowtie2_Index/GRCh38.bowtie2 -I KPGP_00001.bowtie2.sorted.markdup.BQSR.bam -O KPGP_00001.HC.vcf.gz

查找出3,413,019个SNP

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

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