【好久没记录了~最近在准备复试的时候发现了远古时期自己做的大作业，虽然漏洞百出（那时候居然会把个体叫做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 基因组！！

#下载了KPGP_00001这位志愿者的全基因组，存于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 allowed (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

④SNP calling

一、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和参考序列要用统一的！！！

有关全基因组测序的文献复现相关推荐

学习全基因组测序数据分析2：FASTA和FASTQ
本文转载自微信公众号解螺旋的矿工,作者为黄树嘉,已获得授权.黄树嘉写了WGS系列的文章,堪称教科书级别的生物信息学习材料.虽然本平台只关注宏基因组领域,但此系列文章知识体系完善.干货满满,是值得每位专 ...
学习全基因组测序数据分析1：测序技术
本文转载自微信公众号解螺旋的矿工,作者为黄树嘉,已获得授权.黄树嘉写了WGS系列的文章,堪称教科书级别的生物信息学习材料.虽然本平台只关注宏基因组领域,但此系列文章知识体系完善.干货满满,是值得每位专 ...
大熊猫源致病大肠杆菌CCHTP全基因组测序及耐药和毒力基因分析
大熊猫源致病大肠杆菌CCHTP全基因组测序及耐药和毒力基因分析邓雯文1,李才武1,赵思越2,李仁贵1,何永果1,吴代福1,杨盛智2,黄炎1,张和民1,邹立扣2 1. 中国大熊猫保护研究中心,大熊猫国 ...
cfDNA（circulating cell free DNA）全基因组测序
参考资料: [cfDNA专题]cell-free DNA在非肿瘤疾病中的临床价值(好) ctDNA, cfDNA和CTCs有什么区别吗? cfDNA你懂多少? 新发现 | 基因是否表达,做个cfDNA ...
【3月30日直播】新冠病毒全基因组测序——Midnight试剂盒及整体解决方案
识别上方二维码或点击「阅读原文」免费报名参加新冠疫情肆虐全球,基于Nanopore测序技术和数据分析在全球感染性疾病防控中的优势充分显现出来.该平台使用灵活.操作简便.产出快速.分析实时等特征为 ...
《全基因组测序WGS数据分析——1.DNA测序技术》学习笔记
WGS(Whole Genome Sequencing) 指将物种细胞里面完整的基因组序列全部DNA,检测并排列,此技术几乎能够鉴定出基因组上任何类型的突变. 对于人类来说,全基因组测序的价值是极大的 ...
全基因组测序从头测序(de novo sequencing) 重测序(re-sequencing)
全基因组测序全基因组测序分为从头测序(de novo sequencing)和重测序(re-sequencing). 从头测序(de novo)不需要任何参考基因组信息即可对某个物种的基因组进行测序 ...
为什么 Illumina 最新测序仪能将全基因组测序价格降至 1000 美元？
为什么 Illumina 最新测序仪能将全基因组测序价格降至 1000 美元? via Illumina 最新的 HiSeq X Ten 让 1,000 美元一次的 DNA 测序成为现实下面那篇文章 ...
HiFi全基因组测序技术与实例｜HiFi基因组组装软件推荐
HIFI技术的简介 HiFi reads(High fidelity reads) 是Sequel II 三代测序平台推出的兼顾长读长和高准确度的测序序列,一般采用CCS(Circular Conse ...

有关全基因组测序的文献复现