018-基因组的那些事儿

Jul 30, 2018 32 min read cnposts

刘小泽写于18.7.30 偶然间翻出来了18年学习jimmy的”直播我的基因组系列“所做的一些理解，文章写于18.7.30，因为当时感觉工程浩大，所以迟迟没有发出来，但现在我想，“攒着攒着就烂了”，好的内容不能浪费，不踏出第一步，之前的努力就都白费了。当然自己可能还有一些理解不到位的地方，后续内容会补充完整

关于基因组

正常人都是有22对常染色体加XY性染色体。基因组是指生物体所携带的一套完整的单倍体序列，也就是22条+X+Y。每个染色体包括全套基因和间隔序列。他们由A、T、C、G碱基组成，总共长度大约是30亿个碱基。

关于基因检测

随着社会的发展，人们对于健康愈发重视，开始涌现了大量的基因检测，它的个性化定制再加上后续的医师指导，更加准确和便捷获得自身健康信息，预计未来会代替传统体检。

基因检测是在分子水平上对人体遗传密码进行破译，通过单核苷酸多态性和GWAS的分析对人体患病风险进行预测，从而进行预防干预及个体化治疗。目前全基因组测序成本（30X）已经不足一万元，这种测序就是来检测全部的30亿个碱基对是如何排列的，得到从第一个到第30亿个碱基的排列方式。

全基因组检测帮助确诊引起某个疾病的病因，尤其是癌症病人；或者指导有家族性后发遗传病的病人进行有针对性的治疗，比如安吉丽娜·朱莉接受预防性的双侧乳腺切除。

怎么测： 最常用illumina的二代测序，测序长度在150-250bp，取几百万的细胞破碎后，把所有的染色体随机打断成小片段，一个个进行测序，会测得上亿个片段

【还有一种是三代测序，不需要PCR过程，直接对每一条DNA分子进行测序，长度1w-5w nt（因为没有经过PCR，一直是单链状态测，所以不存在碱基对bp，只能称之为碱基nt），准确度要低一些】

测哪里： 也就是测序的样本从哪里获得？

唾液？：唾液肯定可以提取出DNA，而且也最方便。但是会混在口腔微生物的DNA，即使后来通过比对人类参考基因组来去除污染，但最后大概三成数据是要被浪费的。目前基于取唾液兴起的基因检测是测一部分高频变异位点，那不是做的全基因组测序，是利用基因芯片技术进行，成本在三位数

23魔方

血液？除非提供者正患有菌血症（外界的细菌经由体表的入口或是感染的入口进入血液系统后，在人体血液内繁殖并随血流在全身播散），一般血液是最纯净的。从血液里面分离白细胞然后提取DNA的技术也是非常成熟的。

测序报告：

处理流程

数据来源：

一般推荐：全基因组测序，覆盖度30X，也就是90G的raw data，测序策略是PE150，采用illumina的HiSeq X，DNA小片段文库(350bp)进行建库。

几个名词：

覆盖度30X：平均下来能把身体内的30亿个碱基每个都测到30次，因为测序是随机的，必然有一些测序深度高一点，有些低一点
这个30的标准怎么定的？为什么不是20X或者更高的40、50X？有研究做过饱和度分析～看看5~60X的模拟梯度对寻找遗传变异的能力差异大小，结果发现平均深度达到30X的时候，可以覆盖基因组的95%；另外测序深度越高，价格越贵，30X的高性价比足够挖掘到一定量的遗传变异Sequencing depth and coverage: key considerations in genomic analyses - Nature Reviews (2014)
90G raw data：测序深度30X，人类基因组大约30亿碱基，而一亿10^8^ 就等于1Gb的测序数据；拿到的就是3Gb*30X=90Gb。【注意这里的Gb是测序字符的数量】
测序策略PE150:也就是标准的双端测序模式（Paired End），目前双端比单端价格还要便宜，而且一条序列这边测一次，另一边测一次，更准确。所以一般分析基本也没有用单端的了。150就是这边测150bp，那边测150bp。【当然打断的片段一般是大于300bp的，所以每个reads中间会有一部分测不到，这就对了！毕竟reads是随机打断，也就是打断的位置不同。虽然这一条reads的中间部分区域测不到，但是另外的reads就能测到。如果说，一条reads长度200却采用双端150bp，那么中间就会有重叠区域，被测了两遍，这在高通量测序中是非常浪费资源的，每次测都是要花钱的啊！】
Hiseq ：美国Illumina公司作为二代测序仪生产领先企业，自2006年进军基因测序市场以来，陆续发布了HiSeq，MiSeq，NextSeq，NovaSeq等一系列测序仪器。

Hiseq系列～HiSeq 2000，HiSeq 2500，HiSeq 3000，HiSeq 4000 HiSeq系列测序仪问世以来，以通量高，产量大，生产规模著称，能够快速、经济的进行大规模平行测序，在大型全基因组测序，全转录组，全外显子组测序，靶向基因测序方面优势明显。HiSeq 3000/4000系统基于成熟的HiSeq 2500系统，采用创新的有序流动槽技术最大限度提高效率，3.5天内可完成12个基因组、100个转录组或180个外显子组测序

HiSeq X系列——HiSeq X Five，HiSeq X Ten HiSeq X Ten系统的问世完成了人类历史上一大里程碑事件——千元基因组时代的到来。HiSeq X Ten系统是由一套共10台超高通量的HiSeq X仪器组成，其中每台仪器可在3天内产生高达1.8 Tb测序数据，即每天高达600 Gb。10台联合工作，每年能带来超过18,000个人类基因组，而每个基因组的价格约为1000美元，让癌症和复杂疾病的研究达到新的水平

至于NovaSeq嘛，应该是17年开始交付使用，被称为“史上最贵洗衣机”的NovaSeq6000，以其酷炫的外形和美丽的价格（100w美金）成为了高端测序领导者，旨在冲刺“100美元基因组测序”。它的通量更高，运行周期48小时，2个flowcell每次产生大于2Tb的数据。另外还有它兄弟Novaseq 5000，差异就是他们的流动槽，5000可以运行S1、S2两种，6000可以运行S1、S2、S3、S4四种，一个S4流动槽每次运行可达到80-100亿数量的reads / clusters。双S4流动槽运行可以不到两天内解码48个人类基因组（6万亿硷基通量），比双S2流动槽通量翻三倍

DNA小片段（350bp）建库：根据公司不同，将DNA用超声波随机打断成一定长度（如350bp），加接头，作为测序前的准备工作，
Gb与GB你混了吗：
Gb是测序中的数据量，1 Gigabase= 十亿碱基。人类全基因组测序得到了90G的原始数据，也就是900亿碱基。原始数据是fastq格式，而fastq格式是这样的：第二行中一个碱基对应第四行中的一个测序质量
得到的900亿碱基，也对应900亿个质量值，加起来就是1800亿个字符。第一行是测序说明，一般是45个字符，也就是说，每一条测序reads中第一行就有大概45个字符。
那么多少条reads呢？根据PE150计算：测序策略是一条reads包括150bp，现在900亿碱基，就对应900亿/150=60亿条reads 。因此第一行总字符是：60亿*45=270亿个字符。注意到fastq文件共四行，其中1、2、4行的总数量分别为270亿、900亿、900亿，第三行就是一个+，基本可以忽略不计。加起来总共2070亿字符。计算机中，根据编码规则不同，字符与字节对换关系不同。
Fastq文件是ASCII编码文件，其中每一个字符就对应一个ASCII码，也就等于一个字节。计算机的1 GB（Gigabytes） 是1024^3^ 个字节因此，二者对换关系就是：全基因组测序的90Gb对应（2070*10^8^ /1024^3^ ）=193GB计算机存储空间。
或者更快的计算： 测序报告会给出reads数，如果测序策略是PE150，那么占用硬盘空间大小就是n(reads)*(150+150+45)/1024^3^
另外，测序仪下机后的数据都是用gz压缩后的文件.fastq.gz，能压缩2.7倍，大概71G左右。

基本思路：

**科研分析：**测完了基因组，会得到大量的数据，最有价值的分析就是通过与参考基因组的比对，来找出不一样的地方，得到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个要么不准确，要么有点难度】

实际上手

相信你已经磨拳擦掌，想要自己分析一波原始的人类基因组数据，并自己挑出变异基因，这感觉很像一个侦探，逐步深入案情，十分吸引人！

1 计算资源准备

最低配置Linux 8线程+16G+1T硬盘【因为原始数据很大，参考的各个数据库也很大】

查看linux 下计算机配置： CPU：cat /proc/cpuinfo | grep "process" | wc -l 内存：free -g 硬盘：df -h

2 原始数据下载(放在~/project/kpgp_wgs/raw_data下)

选取韩国人Korean Personal Genome Project中的一个人的测序数据KPGP-00001

nohup wget -c -r -nd -np -k -L ftp://ftp.kobic.re.kr/pub/KPGP/2015_release_candidate/WGS/KPGP-00001

关于wget这几个参数：-c 是断点续传；-r是递归下载，因为这里提供的下载地址相当于一个文件夹，下面有很多数据；-nd （no-directories ）不创建目录；-np（no-parent）不要追溯到上一级目录；-k（ –convert-links）转换非相对链接为相对链接【相当于linux中的绝对路径，能够直接访问下载的；因为网页中为了让链接容易记住，设置了大量的相对连接，而我们下载相对链接是访问不进去的】；-L（-L, –relative ）仅仅跟踪相对链接
这些参数设置也不是一成不变的，一般需要设置断点-c、递归-r、-nd

3 参考数据下载

所有的下载之前都要新建好明确的文件夹，比如用户是在home目录下的，一开始空空如也。新建几个目录mkdir biosoft(存放软件)、mkdir project（存放项目）、mkdir reference（存放参考数据）

3.1 参考基因组（放在refernce中）

cd ~/reference
# human genome 19（hg19）
mkdir -p genome/hg19  && cd genome/hg19
for i in $(seq 1 22) X Y M; #指定染色体号下载
do echo $i;
# 一般下载UCSC的数据
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr${i}.fa.gz;
done
gunzip *.gz
for i in $(seq 1 22) X Y M; # 指定染色体号拼成整体
do cat chr${i}.fa >> hg19.fasta
done
rm -fr chr*.fa
# hg38一样的道理
# 参考基因组大小压缩包是900M左右，解压完是3G

小疑点，大问题～基因组版本
人类基因组主要存放在三个数据库NCBI、Ensembl、UCSC NCBI：ftp://ftp.ncbi.nlm.nih.gov/genomes/Homo_sapiens/ARCHIVE/
Ensembl：ftp://ftp.ensembl.org/pub/【然后找release版本，最新是release-93】 UCSC：http://hgdownload.cse.ucsc.edu/goldenPath/hg19 【UCSC的hg系列，是目前使用频率最高的基因组。最新是hg38】

3.2 下载好的基因组需要构建索引（放在reference中）

练习使用bowtie2,hisat2和bwa这3个主流比对软件

软件安装先一律使用conda install -p /YOUR/PATH/ soft_name -y -p制定安装路径；-y允许安装

3.2.1 bowtie2

cd ~/reference
mkdir -p index/bowtie2  && cd index/bowtie2
nohup time bowtie2-build --threads 20 ~/reference/genome/hg19/hg19.fa hg19 1>hg19.bowtie2.log 2>&1 &
#time计算运行时间；--threads设置线程
# 1>hg19.bowtie2.log 新建hg19.bowtie2.log，将运行信息存入；
# 2>&1 2代表错误信息，若有报错，将错误信息覆盖到1的文件中
#hg38参考基因组同理

3.2.2 hisat2

cd ~/reference
mkdir -p index/hisat2   && cd index/hisat2 
nohup time hisat2-build -p 20 ~/reference/genome/hg19/hg19.fa hg19 1>hg19.hisat2.log 2>&1 &
#-p设置线程
#hg38参考基因组同理
#也可以下载现有的：ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/
#下载后的文件夹中还包括make_hg19.sh这个程序，运行可以直接下载hg19基因组

3.2.3 bwa

cd ~/reference
mkdir -p index/bwa   && cd index/bwa
nohup time bwa index -a bwtsw -p hg19 ~/reference/genome/hg19/hg19.fa 1>hg19.bwa_index.log 2>&1 &
#-a指定算法bwtsw；-p指定输出名前缀
#hg38参考基因组同理

3.3 下载基因组注释文件

下载之前先建好存放数据文件夹，先下载gencode数据库中的注释文件。hg38修改了hg19的一些错误，但是在提取指定区域时容易出现一个基因对应两个“分身”位置【见参考4】

这里就有一个问题：GTF和GFF一样吗？
GTF是在GFF的基础上发展来的，二者都是\t分隔的9列文件；GFF包含的信息更多，可以有染色体、基因、转录本信息，而GTF(Gene transfer format)主要描述基因和转录本

mkdir -p ~/reference/gtf/gencode && cd ~/reference/gtf/gencode
## GRCh38 https://www.gencodegenes.org/releases/current.html
mkdir GRCh38_hg38 && cd GRCh38_hg38
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gtf.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.2wayconspseudos.gtf.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.long_noncoding_RNAs.gtf.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.polyAs.gtf.gz
## GRCh37 
mkdir ~/reference/gtf/gencode/GRCh37_hg19 && cd ~/reference/gtf/gencode/GRCh37_hg19
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/GRCh37_mapping/gencode.v28lift37.annotation.gtf.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/GRCh37_mapping/gencode.v28lift37.metadata.HGNC.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/GRCh37_mapping/gencode.v28lift37.metadata.EntrezGene.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/GRCh37_mapping/gencode.v28lift37.metadata.RefSeq.gz

然后按照我们的需要提取出一些指定区域，比如蛋白质编码区，长链非编码RNA区、假基因区(pseudogene是与正常基因相似，但丧失正常功能的DNA序列)以及hg19、hg38各自的全基因位点区域信息

## 提取GRCh37(hg19)全基因信息
cd ~/reference/gtf/gencode/GRCh37_hg19
zcat gencode.v28lift37.annotation.gtf.gz  | perl -alne '{next unless $F[2] eq "gene" ;/gene_name \"(.*?)\";/; print "$F[0]\t$F[3]\t$F[4]\t$1"}'> allgene.hg19.position

对这一行代码的解读：zcat是在不解压缩文件的前提下查看文件，结果通过管道符|传给后面命令；perl接受了gtf文件中一行一行的信息，要统计各个基因所在染色体、起始位置、终止位置、基因名这四种信息，就要先对gtf文件做解读
处理tab分割的文件的每一行内容，perl -alne是必备佳品：例如：while(<>){chomp;@F=split /s+/,$_;print "$F[2]\n"}也就相当于 perl -alne 'print $F[2]～多么简洁！ -e打开单行输入模式；-n打开循环模式；-l对输入内容自动chomp，对输出内容自动添加换行符；-a自动分隔模式 那么这里的{next unless $F[2] eq "gene 意思就是如果第二列（Feature信息）不等于“gene”就看下一行； /gene_name \"(.*?)\";/ 是将gene_name后面引号中(注意引号要用转义一下)的部分存入待引用变量中，那么存入多少呢？这里的(.*?)就起作用了：.是除了\n的任意字符，*是取之前字符的0个或者n个。perl默认使用贪婪模式匹配，也就是如果不限制，他能一直匹配到结尾，因此?就是限制贪婪模式的，这三个连用，就命令perl，只准匹配到后一个引号，不能继续了哈！接下来就是打印输出了，print "$F[0]\t$F[3]\t$F[4]\t$1" 将第一列、第四列、第五列再加上前面待引用部分【用$1调取】打印出来，中间用分隔符分隔

要知道提取什么位置，就必须要先了解gtf的内容，除了之前说的那几列，有一个很重要的知识点：gene type的类型重复多次出现，有什么规律吗？

##统计hg19编码蛋白的基因信息
zcat gencode.v28lift37.annotation.gtf.gz | grep protein_coding | perl -alne '{next unless $F[2] eq "gene"; /gene_name \"(.*?)\";/;print "$F[0]\t$F[3]\t$F[4]\t$1"}'>hg19_protein_coding.position

##统计hg38的长链非编码RNA信息
zcat gencode.v28.long_noncoding_RNAs.gtf.gz| perl -alne '{next unless $F[2] eq "gene"; /gene_name \"(.*?)\";/;print "$F[0]\t$F[3]\t$F[4]\t$1"}'>hg38_lncRNA.position
##统计hg38的假DNA信息(第三列type信息全是transcript)
zcat gencode.v28.2wayconspseudos.gtf.gz | perl -alne '{next unless $F[2] eq "transcript"; /gene_name \"(.*?)\";/;print "$F[0]\t$F[3]\t$F[4]\t$1"}'> hg38_pseudos.position
##统计hg38的全基因信息
zcat gencode.v28.annotation.gtf.gz | perl -alne '{next unless $F[2] eq "gene"; /gene_name \"(.*?)\";/;print "$F[0]\t$F[3]\t$F[4]\t$1"}'>allgene.hg38.position

结果得到了hg19的编码蛋白基因信息

要查找基因信息只需要grep基因名就好，标准基因名可在https://www.genecards.org/中查找，例如要查找抑癌基因TP53和乳腺癌基因BRCA1的信息，就可以

grep TP53 hg19_protein_coding.position
grep BRCA1 hg19_protein_coding.position

4 质量控制(放在kpgp_wgs/qc下)

使用fastp进行质控，速度十分快，我之前测试了24G的双端测序文件，不到六分钟处理完成。下载的KPGP原始数据共68G，6个双端测序文件，平均一个用时1300秒

因为是从数据库中下载的现成的测序文件，质量已经调到了最好，用fastp质控后，改动不是很大【如果是测序质量差的文件，在fastp过滤后会有明显差别】如图是我们的KPGP的原始数据

可以自己统计GC、Q20、Q30

#!/usr/bin/perl
opendir (DIR, "./") or die "can't open the directory!";
@dir = readdir DIR;
foreach $file ( sort  @dir) 
{
next unless -d $file;
next if $file eq '.';
next if $file eq '..';
$total_reads=  `grep '^Total' ./$file/fastqc_data.txt`;
$total_reads=(split(/\s+/,$total_reads))[2];
$GC= `grep '%GC' ./$file/fastqc_data.txt`;
$GC=(split(/\s+/,$GC))[1];
chomp $GC;
open FH , "<./$file/fastqc_data.txt";
while (<FH>)
	{
	next unless /#Quality/;
	while (<FH>)
		{
		@F=split;
		$hash{$F[0]}=$F[1];
		last if />>END_MODULE/;
		}
	}
$all=0;$Q20=0;$Q30=0;
$all+=$hash{$_} foreach keys %hash;
$Q20+=$hash{$_} foreach 0..20;
$Q30+=$hash{$_} foreach 0..30;
$Q20=1-$Q20/$all;
$Q30=1-$Q30/$all;
print "$file\t$total_reads\t$GC\t$Q20\t$Q30\n";
#print "$all\n";
# 
}

5 比对：终其一生，只为遇见你(放在kpgp_wgs/align下)

何为终其一生？ 每条测序reads长度是150bp，只有reads这些ATCG字母，对我们是没有用的，起码要知道这些ATCG的排列在基因组的什么地方吧。人类的参考基因组30亿个碱基，如果要150一段150一段去观察的话，真的，一条reads要终其一生，只为遇见那个和他匹配的那个她。更何况，一个全基因组测序下来就有至少6亿条reads😱

但！我们有不怕苦不怕累的计算机！可以使用软件bwa或者bowtie，将fastq测序文件与参考基因组做对比，得到的结果是sam格式的文件【如下图】，其中可以看到每条reads在参考基因组的位置，在哪一条染色体，或者在这条染色体的哪个位置

5.1 Fastq to sam

为了得到这个结果，我们可以使用bwa软件进行操作

#bwa mem -t线程数 -M处理同一个reads比对到参考基因组上不同位置的情况 bwa_index（要定位到前缀，不能只定位到文件夹名，比如这里倒数第二个hg19是文件夹名，最后一个hg19才是前缀）输入fastq1 输入fastq2 1>标准输出（sam格式）2>标准错误输出
for i in $(seq 1 6);do
bwa mem -t 10 -M ~/reference/index/bwa/hg19/hg19 ../raw_data/KPGP-00001_L${i}_R1.fq.gz ../raw_data/KPGP-00001_L${i}_R2.fq.gz 1>KPGP-00001_L${i}.sam 2>KPGP-00001_L${i}.bwa.align.log
done

【使用60线程运行，平均50分钟比对完一个】得到的结果是sam格式的，sam的全称是sequence alignment/map format，往往文件很大，为了节省空间又保证同样效果，他的孪生兄弟bam诞生【bam是sam的二进制文件，占用空间要小很多】

5.2 sam to bam

需要用samtools进行格式转换【转换依旧使用60线程，大概半个小时就处理好了，可以看出bam比sam节省了三分之二的空间】

#sam2bam -b设置输出为bam格式 -S 指定输入为sam
samtools view -b -S file.sam > file.bam
#bam2sam -h意思是加上header
samtools view -h file.bam > file.sam

5.3 sam/bam文件格式

有了文件，以后的大部分分析都是基于sam格式进行的，因此理解格式至关重要

第3和第7列，可以用来判断某条reads是否成功比对到了基因组的染色体上，另外两条PE reads是否比对到同一条染色体；第1，10，11列可以还原成测序数据的fastq格式

5.4 对bam文件进行排序=>sorted_bam

for i in `seq 1 6`;do
samtools sort  -@ 15 -o KPGP-00001_L${i}.sorted.bam KPGP-00001_L${i}.bam
done

每个bam文件按照最左侧比对位置的先后进行排序，方便接下来的构建索引

5.5 对sorted_bam文件构建索引=>sorted_bam.bai

for i in `seq 1 6`;do
samtools index KPGP-00001_L${i}.sorted.bam
done

bam文件进行默认情况下的坐标排序后,才能进行index。目的是为了更快地访问bam文件，结果会生成.bai格式的索引文件

samtools index是针对bam用的，如果要对sam文件构建索引，要用samtools tabix命令

5.6 查看参考基因组每个位点的比对情况=>mpileup.txt

利用samtools的mpileup功能可以生成bcf文件，之后可以结合bcftools进行SNP和InDel的分析；

参数：-f：输入有索引的参考基因组fasta；

-g：输出二进制的bcf格式（不加-g就只生成文本文件，统计了针对不同的测序比对文件，参考序列中每个碱基位点的情况，每一行表示参考序列中某个碱基的比对结果）

for i in `seq 1 6`;do
samtools mpileup -f ~/reference/genome/hg19/hg19.fa KPGP-00001_L${i}.sorted.bam > KPGP-00001_L${i}.mpileup.txt
done

查看各个测序样本bam比对结果

for i in `seq 1 6`;do
samtools flagstat KPGP-00001_L${i}.sorted.bam >KPGP-00001_L${i}.flagstat.txt
done

以上几步操作每个测序比对文件会形成

5.7 将各个bam文件融合=>merge.bam

原来只是各个测序文件比对的结果，他们虽然都sort和index了，也有了各自与参考基因组的位点比对情况。但是我们不能从整体角度去查看，要想统一原本的“各自为政”，我们需要进行merge

samtools merge -@ 20 KPGP-00001.merge.bam *.sorted.bam

5.8 对merge后的bam也进行排序=>sorted.merge.bam

##sort
samtools sort -@ 20 -O bam -o KPGP-00001.sorted.merge.bam  KPGP-00001.merge.bam
##index
samtools index KPGP-00001.sorted.merge.bam

5.9 对merge_bam进行拆分=>n个小bam

结果得到的KPGP-00001.sorted.merge.bam是61G，这样的文件是不可能直接导入IGV查看的，有没有什么方法可以缩小文件进行查看呢？

查看一个大文件不靠谱，那么可不可以将大文件拆分开来，拆成一个个的小文件呢？利用bamtools split是可以的。另外，如果程序随意给我们排序好的bam文件给拆了，那么之前的排序工作也没有任何意义，因此，我们可以指定参数-reference，让程序根据参考序列（也就是根据bam第三列的染色体编号） 进行拆分

bamtools split -in KPGP-00001.sorted.merge.bam -reference

下面这张图就是合并bam、合并后排序的bam以及拆分后的bam全家福

5.9.1 【比对不上的】统计未比对上的测序reads=>unmapped.counts

**方案一：**上面我们得到了一个文件KPGP-00001.sorted.merge.REF_unmapped.bam ，统计具体比对信息使用一个小脚本

samtools view KPGP-00001.sorted.merge.REF_unmapped.bam | cut -f 3,6 |sort -V|uniq -c >unmapped.counts
## 将bam文件第3列、第6列提取出来,也就是参考序列这一列，还有CIGAR列
## 【CIGAR是*，才能说明左右端是匹配失败的】
## sort将提取出来的内容排序，-V表示按照正常阅读方式显示（默认按照ASCII码方式排序，例如，chr1下一个是chr10而不是chr2）
## uniq 将相同的归为一类， -c统计这一类有多少个

但是发现，结果只有一行：1931716 * * 也就是只显示了未比对上的总数，并没有按期望显示各个染色体分别的统计

关于解释：【CIGAR是*，才能说明左右端是匹配失败的】看下图👇

**方案二：**如果想看根据每一条染色体进行的统计，可以用bamtools split 的-mapped 参数，它会将比对上的和未比对上的分开，再将未比对上的bam中第三列refname和第六列CIGAR取出来

bamtools split -in KPGP-00001.sorted.merge.bam -mapped #这样就是按整体划分为比对上和未比对上，一会直接看未比对上的，那里就包括所有的染色体
samtools view KPGP-00001.sorted.merge.UNMAPPED.bam | cut -f 3,6 |sort -V|uniq -c >unmapped.counts

结论： 方案一只显示总的结果；方案二才能显示各个染色体上存在未比对的具体情况

5.9.2 【比对不同的】统计左右两端测序数据比对到不同染色体的reads=>unpaired.counts

在测序过程中，配对的reads来自同一个DNA，一般可以比对到同一个染色体
sam 文件中第三列和第七列表示了read1、read2比对到了哪条染色体。如果第七列为“=”，那么两条read比对到了同一个染色体，可以认为它们是配对的；
但往往reads能比对到不同染色体上，会更有意义，揭示了融合基因的可能性或者说本来有两个基因就很相似【融合基因指两个基因的全部或一部分的序列相互融合为一个新的基因的过程，可能是染色体易位、中间缺失或染色体倒置所致】
因此，如果要找比对到不同染色体的reads，只需要找第七列不是“=”的情况

方案一： 仅对比对上的reads做统计

samtools view KPGP-00001.sorted.merge.MAPPED.bam  | perl -alne '{print if $F[6] ne "=" }'>unpaired_1.sam
cut -f 3,7 unpaired_1.sam | sort | uniq -c > unpaired_1.txt

**方案二：**对全部的reads做统计

samtools view KPGP-00001.sorted.merge.bam  | perl -alne '{print if $F[6] ne "=" }'>unpaired_2.sam
cut -f 3,7 unpaired_2.sam | sort | uniq -c > unpaired_2.txt

结论： 方案二比一多了【两条read都没比对上】的情况

5.9.3 【多重比对的】找出multiple-mapping

什么是multiple-mapping（多重比对）：如果一条或一对 read在基因上有多个比对位置，BWA从中选出一条最好的；如果有多个一样好的比对位置，随机选取一个。
后期检测SNP或InDel的时候需要高质量的比对结果，在sam文件的第五列是MAPQ值，表示比对质量，数值越大质量越高。在IGV中，每一个矩形箭头都是一条比对好的reads，白色的箭头说明比对质量值为0，灰色的是正常的reads。
BWA软件处理策略与其他比对软件不同，它认为大于1的就是unique mapping。如果要提取multiple mapping，那么提取MAPQ值为0即可。

samtools view KPGP-00001.sorted.merge.bam| perl -alne '{print if $F[4] eq 0}' > multiple_mapping.sam
##结果得到的文件大小为14G

5.9.4 分别统计覆盖度和平均测序深度

首先理解什么是覆盖度和测序深度？ 覆盖度是匹配到染色体上reads的数目/染色体长度【先统计sort后的bam文件中，有多少reads比对到了chr1，有多少在chr2…，把这些数字分别和对应的染色体长度相除，得到一个百分数。解释一下就是：有百分之多少的reads是能比对到哪条染色体的】
测序深度 是染色体每个位点得到的reads总数/染色体长度【先统计sort后的bam文件，结果中包括了每条染色体的各个位点比对到了多少条reads，也就说明这个位点被测得了多少次。举个🌰，下面图中第一个“野生的命令”中，红框框住的就是chr1上的9999位点发现了17条测序reads，10000位点发现了18条测序reads…我们按照不同染色体将每一个位点比对到的reads加起来，就是分子；分母和覆盖度一样，就是染色体长度。如果算出来的chr1结果是33，那么就说明chr1的测序深度是33X】平均测序深度就是把所有染色体的测序深度求平均值
主要是用samtools的depth和mpileup功能将bam文件转为统计数据，然后再用脚本进行分别统计～

先比较几个bam转为统计数据的处理方案：

上面的介绍，知道了计算覆盖度和平均测序深度的分母都是各个染色体长度，下面就用一个perl脚本分别计算分子

nohup samtools mpileup KPGP-00001.sorted.merge.bam| perl -alne '{$pos{$F[0]}++;$depth{$F[0]}+=$F[3]} END {print "$_\t$pos{$_}\t$depth{$_}"foreach sort{$a <=> $b}keys%pos}' 1>coverage_depth.txt  2>coverage_depth.log &

解读脚本：nohup…&是将程序放后台进行，防止中途断掉； samtools mpileup KPGP-00001.sorted.merge.bam是将sort后的bam文件转化为统计数据，第一列是reference name，这里就是染色体号，第四列是每个位点比对到了多少条reads； perl -alne '{}' 是处理tab格式文件常用格式，$F[0]就是第一列的染色体号，{$pos{$F[0]}++ 统计有多少chr1，多少chr2…，并储存在数组变量pos中，$depth{$F[0]}+=$F[3]将第四列的结果累加到数组变量depth中； END表示前面计算结束，开始打印{print "$_\t$pos{$_}\t$depth{$_} 。每运行一次，都会产生三个值：$F[0]、{$F[0]}++、$F[0]}+=$F[3]储存在内存中，于是可以通过$_调用，定义的输出格式为染色体号\tpos值\tdepth值，中间使用tab分隔； foreach定义循环结构，对每一行进行排序sort{$a <=> $b}keys%pos} ，其中$a <=> $b 是飞碟符号，表示从小到大排序，然后按照数组pos中的keys值也就是chr的值进行排序；【补充：如果要逆序，用$b <=> $a】最后结果标准输出到coverage_depth.txt，错误信息到coverage_depth.log

16_2

可以看到，Y染色体的覆盖度和测序深度少的离谱，21、22染色体的覆盖度也比较低。原因是什么呢？是数据质量真的这么差么？

当然，也可以换一种转换bam为统计数据的工具，例如可以再用samtools depth再统计一遍，测序深度应该会有提高，但也不是很理想。

这个原因就是： 我们使用的第三列数据ch_len，是利用samtools view -H KPGP-00001.sorted.merge.bam 产生的表头信息。而这个bam又是和参考基因组比对得到的。因此，我们统计的时候并没有对参考基因组进行任何处理，除了ATCG碱基之外，还有一类存在，那就是N ，计算覆盖度的时候，我们把大量的N加入了分母中，造成数据偏低

重新整理数据 ：统计N的个数，然后将N从总数中去掉

##如果使用的参考基因组是human_g1k_v37.fasta，根据它的格式>后面是数字，没有chr三个字符,比如 >1 dna:chromosome 
less human_g1k_v37.fasta | perl -alne '{if (/^>(.*?)\s/){$chr=$1}else{$N_count{$chr}+=($_=~tr/N//)}}END{print "$_\t$N_count{$_}"foreach sort {$a <=> $b}keys %N_count}' > hg_v37_chr_len.txt
#$1指的就是小括号中的匹配内容
#$_ 是undef,相当于一个随时待命的，这里就是$chr
#根据v37的格式，大于号开头，接一个数字，然后空格后又有具体信息
#我们想要数字，用.*，并且加上?防止贪婪匹配，还要指明匹配的内容在空格(\s)之前
#foreach sort {$a <=> $b}按chr的数字从小到大排序，排序的主体就是keys %N_count，keys是键=》$chr，%N_count是值=》$N_count

##如果使用的是参考基因组hg19/hg38.fasta，他们的格式是>后是chr加数字，比如>chr1
less hg19.fa | perl -alne '{if (/^>(.*)/){$chr=$1}else{$N_count{$chr}+=($_=~tr/N//)}}END{print "$_\t$N_count{$_}"foreach sort {$a <=> $b}keys %N_count}'

结果好了很多，平均测序深度基本能到30X，但是Y确实比较低

16_3

16_4

5.9.5 作出详细的报告列表=>reads+不同depth下的coverage信息

使用samtools flagstat功能统计reads的一些信息

##用merged_bam做一个总体的
samtools flagstat -@ 20 KPGP-00001.sorted.merge.bam >KPGP-00001.sorted.merge.flagstat.txt

使用mileup再加上perl小脚本统计不同测序深度时的覆盖度信息 目的：得到全基因组区域中碱基覆盖深度不低于多少X的比例，如：Coverage_at_least_4X、Coverage_at_least_10X、Coverage_at_least_20X并且能做出来测序深度与覆盖度关系图

samtools mpileup KPGP-00001.sorted.merge.bam | perl -alne '{if ($F[3]>100){$depth{"100 plus"}++}else{$depth{$F[3]}++}}END{print "$_\t$depth{$_}"foreach sort{$a <=> $b}keys%depth }' > all_id_depth.txt

和5.9.4的脚本很类似，也是转换数据-》计算-〉打印-》排序。不同之处就在于perl的计算过程，以下解读脚本：

if ($F[3]>100){$depth{"100plus"}++} 如果第四列大于100，也就是指某个位点比对的reads数超过一百（即该位点测序深度大于100），那么就将累加的每个位点测序深度频数赋给“100 plus”；如果小于100，单独统计每个测序深度的比对reads数

作出测序深度与测序位点数的关系图

# 测序深度与测序位点数的关系 作图代码
a=read.table("all_id_depth.txt",stringsAsFactors = T)
ggplot(a, aes(x=1:nrow(a),y=V2,col="orange"))+geom_line(size=1.5,alpha=0.5)+
  labs(x="Depth",y="Positions")+geom_vline(xintercept = 30,color="purple",size=0.4)+
  scale_x_continuous(limits=c(0,100),breaks=seq(0,100,5))+
  scale_y_continuous(limits=c(80000,130000000),breaks=seq(80000,130000000,40000000))+
  theme(axis.text=element_text(size=16),axis.title=element_text(size=20,face="bold"))

作出不同测序深度的累计分布，也就是看有多少大于10X、20X、30X

# 不同测序深度的累计分布 作图代码
a=read.table('all_id_depth.txt',stringsAsFactors = F)
over100=a[1,]
a=a[-1,]
a=rbind(a,over100)
a[,3]= cumsum(as.numeric(a[,2]))
a[,4]= 1-a[,3]/sum(as.numeric(a[,2]))
ggplot(a,aes(x=1:nrow(a),y=V4,col="orange"))+geom_line(size=1.5,alpha=0.5)+
  labs(x="Depth",y="Percentage")+
  geom_vline(xintercept = 10,color="purple",size=0.4,linetype="dashed")+
  geom_vline(xintercept = 20,color="navy",size=0.4,linetype="dashed")+
  geom_vline(xintercept = 30,color="red",size=0.4,linetype="dashed")+
  scale_x_continuous(limits=c(0,100),breaks=seq(0,100,5))+
 theme(axis.text=element_text(size=13),axis.title=element_text(size=15,face="bold"))

当然这是最整个基因组进行的统计，如果要统计每一条染色体的测序深度信息，需要一个循环脚本，然后统计结果直接将上面代码中txt替换掉就可以

ls KPGP-00001.sorted.merge.REF*.bam | while read file
do
#echo $file 测试一下能否输出全部的小bam的名字
samtools mpileup $file | perl -alne '{if ($F[3]>100){$depth{"100 plus"}++}else{$depth{$F[3]}++}}END{print "$_\t$depth{$_}"foreach sort{$a<=>$b}keys%depth}' > $file_depth.txt
done

这里被困住了两天，需要学脚本，回忆R和ggplot2的知识，终于算是搞定了，继续出发！然后最近办小型培训班也需要付出精力，有点累但值得！写于18.8.9

前面都是全基因组水平或者个染色体水平的统计信息，主要和测序质量有关，接下来开始分析真正的变异信息

6 练手-初步检查几个基因变异位点

搜索人体内有趣的基因，就会出现许多好玩的基因，让我们看一下我们研究的样本有没有这样的基因吧，可以由此推断一个人的特性哦，就像侦探一样，有趣！

#上网搜索一下“人体内有趣的基因”就会有许多结果
V1aR基因是雄性标志性出轨基因; =>标准命名“AVPR1A” 
FOXO3A长寿基因； =>标准命名“ FOXO3”
5-HTTLPR快乐基因；=> 标准命名“ SLC6A4”
hDEC2不眠基因 => 标准命名“ BHLHE41”

首先使用GeneCards数据库查找标准的基因命名 ，比如V1aR基因
然后将各个基因坐标信息添加到一个文本文件中

从merge的全基因组bam文件中根据感兴趣基因列表，提取他们的变异情况

cat list_genes|while read id;do
chr=$(echo $id |cut -d" " -f 1 )
start=$(echo $id |cut -d" " -f 2 )
end=$(echo $id |cut -d" " -f 3 )
gene=$(echo $id |cut -d" " -f 4 )
echo $chr:$start-$end  $gene
samtools mpileup -r $chr:$start-$end -ugf ~/reference/genome/hg19/hg19.fasta ./bamFiles/KPGP-00001.sorted.merge.bam | bcftools call -vmO z -o $gene.vcf.gz
done

参考：

Illumina测序仪https://www.sohu.com/a/143929658_603295
Hiseq Xhttp://www.illumina.com/content/dam/illumina-marketing/documents/products/appnotes/appnote-hiseq-x.pdf
hg19 (GRCh37) 与 hg38 (GRCh38) 数据差异比较https://wenku.baidu.com/view/eaa57f642bf90242a8956bec0975f46527d3a7a4.html
关于一个基因在hg38基因组上有两个定位的问题：http://www.bio-info-trainee.com/1991.html
Perl单行命令参数讲解：https://www.jb51.net/article/123326.htm
Sam/Bam文件Flag值转换：https://broadinstitute.github.io/picard/explain-flags.html
关于samtools：https://www.jianshu.com/p/240bffd7cf00
关于长寿基因FOXO3A：https://singularityhub.com/2010/02/19/want-to-live-forever-better-hope-you-have-the-right-foxo3a-gene/#sm.00015qhufi51pfhxr9u2p7e56y57l
GeneCards数据库：https://www.genecards.org/
基因名标准化数据库（HGNC）：https://beta.genenames.org/

WGS