生物信息前期入门大纲

Posted 2020-11-26 hongwan

零、前言

“不懂就问”，但是这里的“问”是指问百度和谷歌，实在不懂再问人!

上谷歌教程，参考文章：谷歌浏览器插件与电脑软件推荐

(1)生信论坛推荐
生信技能树（生信菜鸟团），有很多优秀的帖子，以及对应的微信公众号，其它生信微信公众号优秀的文章它会转载，所以关注这个够了，历史消息值得翻阅一遍；
PLoB:https://www.plob.org/
Biostars，国外生信论坛

(2)文本文件浏览软件
nodepad++、sublime text、破解版的UE

一、linux操作

Linux 是操作系统的一种，我们台式机笔记本用的系统是Windows，win7，win10等，苹果系统叫mac os,手机操作系统常见的有苹果的ios，android，可以用于运行其它软件；
但是我们使用linux一般使用的不是图形界面，不用鼠标操作，而是使用命令行，不同单词表示不同的指令；比如：ls 表示列出当前目录所有文件，cd 表示切换目录，等等；使用linux主要是因为生信软件很多都是Linux的,而且可以批量操作；

常用Linux指令.sh记下来，遇到不懂的到时候再百度；
推荐书籍：《鸟哥私房菜》；后续当字典用；

(1)如何登陆linux
软件：暂时推荐Xshell；点击：文件-新建-输入主机ip-确定-输入用户名-确定-密码-确定

之后点击：打开-选择主机即可；
字体太小，调大字体：

还可以使用：FillZilla,可视化，主要用于传输文件，把文件上传到服务器或者下载到本地
·登录方法：
选择“文件-站点管理器-协议选择SFTP--输入账号密码登录”

Linux最常见错误：
No such file or directory：找不到文件，没有这个文件或者目录写错了；
其它请百度这个错误；

投递任务时，请注意内存使用； 节约硬盘：原始文件获取到结果之后请删除，或者转为bam文件，删除sam和fastq文件； 保存好每一步代码，因为你很大可能会重跑

二、R语言和Rstudio

R语言是常用于统计分析、绘图的编程语言，入门简单，适合没有编程基础的学习；
书籍途径：《R语言实战》
请先安装R以及Rstudio，Rstudio 是R语言的IDE，更加方便使用R编写代码；
最主要了解：读取文件、调用函数、绘图、保存文件；

Rstudio有四个页面：左上为代码编写框，左下为运行框：显示运行的代码和输出结果以及报错信息等，直接输入回车即运行，右上位变量情况和历史命令；右下包含有图形显示，包管理，帮助页面等；

R包安装：
选择Tools-install packages—输入包名—install

如果没有，比如安装limma：
在bioconductor中搜索，
输入source("https://bioconductor.org/biocLite.R")
,然后输入biocLite（”包名“）

source("https://bioconductor.org/biocLite.R")
biocLite("limma")

遇到新包，请先? / help 查看包的说明页面，并一条条运行其官方例子，可以快速了解此包；

三.python

Python是目前比较火的编程语言，简单易学，常用于文本处理，爬虫和深度学习等处理；
推荐教材：
廖雪峰python教程

可以先看到函数那一部分；
软件安装，请最好安装anocoda,可以少趟坑
编辑器推荐：sublime text

四.测序和芯片数据下载

(1)测序数据

GEO上有四类数据GSM, GSE, GDS, GPL
GSM是单个样本的实验数据
GDS是人工整理好的关于某个话题的GSM的集合，一个GDS中的GSM的平台是一样的
GSE是一个实验项目中的多个芯片实验，可能使用多个平台
GPL是芯片的平台，如Affymetrix， Aglent等

比如：GSM2759564，在ncbi，选中sra,搜索GSM2759564，
可以点击SRR值，点击download，会提示你使用fastq-dump或者Prefetch 下载，但是不推荐，因为老是断；

可以在此页面找到ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/
SRP对应页面的sra文件夹；获取到sra的链接，其实是有规律的，但是不一定有，所有还是需要确认下；
ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP前三位/SRP值/SRR值/SRR值.sra
这里对应的就是：ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP124/SRP124514/SRR5990244/SRR5990244.sra

然后你可以使用wget、axle等常用下载工具下载，或者使用超级快的Aspera下载；
替换链接即可

ascp -QT -l 100M -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh [e-mail]: /sra/sra-instant/reads/ByStudy/sra/SRP/SRP124/SRP124514/SRR5990244/SRR5990244.sra

具体教程：[使用Aspera从EBI或NCBI下载基因组数据]modified(https://blog.csdn.net/likelet/article/details/8226368)

(2)芯片数据下载

Supplementary file页面下方Supplementary file即是；
Series Matrix Files为注释文件，将芯片id与基因名对应起来

五.测序原理

由于生信大部分时间都在于测序数据打交道，了解测序原理更有利于后续分析；
推荐：陈魏学基因，爱奇艺和YouTube都有，很适合入门了解；
先了解mRNA测序，Illumina测序等常用测序的原理

测序技术：从零开始完整学习全基因组测序（WGS）数据分析：第1节 测序技
illumina 测序原理介绍

下文为Illumina测序文字部分：
主要部分：建库——桥式PCR——测序

首先是Flow cell，是测序反应发生的位置，每一个flowcell上都有8条泳道，用于测序反应，可以添加试剂，洗脱等等.
在Flow cell玻璃表面有两种DNA引物序列，是与要测序的DNA文库接头序列互补的。

接着是DNA文库：DNA文库是指许多DNA片段，在两边接上了特定DNA接头。
先通过将基因组DNA打断成一定长度的片段，打断之后会出现末端不平，使用酶补平，并在3‘端使用酶加上一个特异的碱基A，连接酶接上人工特定的接头，方便扩增和测序；然后就进行PCR扩增

然后是核心部分：桥式PCR，就是把文库种到芯片上并进行扩增的过程
构建好的DNA文库 与芯片引物互补杂交；然后加入dNTP和聚合酶，合成新的DNA链，加入NaOH解开双链；不和芯片连接的原来的链被冲走；加入中和液体，DNA一端和第二种引物互补杂交。

再加入酶和dNTP,聚合酶沿着引物合成新链，再加碱解链，再加中和液，两条DNA链都会和新的引物杂交，反复这个过程。DNA呈指数形式增长；形成很多cluster,1个cluster就是一群完全相同的序列

最后一步就是测序：首先通过化学反应，把一个引物的特定基团切断掉，加碱，留下共价键连接的单链；

聚合酶根据互补原理结合,结合对应的dNTP，加入的碱基有两个特点：第一个是脱氧核糖3号位加入了叠氮基团而不是常规的羟基，保证每次只能够在序列上添加1个碱基；另一方面是，碱基部分加入了荧光基团，可以激发出不同的颜色。

一个循环只能延长一个碱基，停止，放到显微镜进行激光扫描，根据发出的荧光判断新合成的碱基是哪个碱基。然后加入化学试剂把叠氮基团和旁边荧光基团切掉，然后进行新一轮的合成。

还有通过index 区分样本来源和双端测序：把原来的链去掉，用测序生成的新链测序，相同的过程。

六、RNA-seq 分析流程

流程跑过一次就基本会了的，只看不跑还是不会的；

主要包含：数据下载、格式转换、数据质控（FastQC、MultiQC）、数据修剪( Trimmomatic 、afterQC)、序列比对(hisat2、TopHat、STAR、Bowtie、Bowtie2)、排序（samtools）
定量（RSEM、htseq）、差异表达基因分析和功能富集分析。

第一步:将sra 转为fastq格式

使用fastq-dump将sra 转为fastq格式文件
fastq格式数据格式主要包含四行：第一行是坐标信息、第二行是序列信息、第三行是附加信息、第四行是序列的质量信息：在测序仪进行测序的时候，会自动根据荧光信号的强弱给出一个参考的测序错误概率（error probility，P），P值肯定是越小越好。将P取log10之后再乘以-10，得到的结果为Q。把这个Q加上33或者64转成一个新的数值，称为Phred，最后把Phred对应的ASCII字符对应到这个碱基，使得1个符号与1个碱基一一对应。

第二步:数据质控

做质控是因为二代测序是边合成便测序的，随着合成链的增长，DNA聚合酶的效率会不断下降，特异性也开始变差，越到后面碱基合成的错误率就会越高;测序数据的质量好坏会影响我们的下游分析。

一般是使用fastqc这个软件，结果部分会包含：总体的基本信息，（平台版本、reads数目、测序长度和GC含量）已经每条序列每个碱基质量、每条序列的测序质量统计、GC 含量统计、序列测序长度统计、序列Adapter、重复短序列（在序列中某些特征的短序列重复出现的次数）；比如第一个代表的是每个序列每个碱基的质量，一般要求大于20。

但是如何样本数目过多：不可能一个个看；可以自己写脚本去统计每一个样本的各个指标；
也可以通过别人的软件；比如MUltiQC，可以统一查看fastqc结果；比如我们的数据，查看数据总体情况：read各位置的N含量；N 含量比较高:意味着，测序的光学信号无法被清晰分辨，测序系统或者测序试剂的错误。还有接头：建测序文库，测序接头就是在这个时候加上的，其目的一方面是为了能够结合到flowcell上，另一方面是当有多个样本同时测序的时候能够利用接头信息进行区分。

详细教程：从零开始完整学习全基因组测序数据分析：第3节 数据质控

第三步:数据修剪

目的就是为了去除数据中质量低的reads和接头序列；
之前是用Trimmomatic,自己通过脚本提取fastqc结果中的adapter，设置参数修剪；
而这次用的是fastp(https://github.com/OpenGene/fastp)，可以自动过滤低质量reads、做序列剪切、消除误差和质量控制；
结果会有Good、bad、和QC 三个文件夹：分别是质量好的、坏的fastq文件和运存结果报告 ;结果文件中有:结果总结：总read和碱基数、过滤的read和碱基数、接头切除数目、过滤之前和之后的read质量曲线、各个碱基含量分布、GC含量。

第四步：序列比对和排序

目的是将无序的read片段 组装回参考基因组，
经过DNA建库和测序之后，FASTQ文件中紧挨着的两条read 没有位置关系，都只是随机来自于原本基因组中某个位置的短序列。因此需要一个个去和参考基因组比较，找到每一条read在参考基因组上的位置，然后按顺序排列好。
所以下载参考基因组和基因组注释文件；
为了提高比对速度，就需要根据参考基因组序列，经过BWT算法（具体看不懂）转换成index；
把read和index进行比较；
然后使用samtools将sam文件转为二进制的bam文件
第一步的比对是按照FASTQ文件的顺序把read逐一定位到参考基因组上的，比对后得到的结果文件中，每一条记录之间位置的先后顺序是乱的，我们后续分析需要按照顺序从小到大排序下来才能进行。所以需要排序。

左边是基因组注释文件数据组成：9列：
右边是常用的比对软件：小片段用bowtie和Bowtie2，小于50bp推荐用1；大于50bp推荐用2；常用的还有Tophat2 ，但是作者推荐HISAT2，更快，比STAR还快；

第五步：Reads计数、合并表达矩阵

目的： 统计read在哪些基因上，以此来计算不同基因的表达量，然而都是得到了相对定量的值，即某一个基因转录本占样本中所有转录本的百分比，所以需要标准化，这里标准化方法有FPKM、RPKM和TPM，目的是为了排除测序深读和基因长度的影响

RPKM：落在这个基因上的总的read数与这个样本的总read数和基因长度的乘积的比值
TPM：先对每个基因的read数用基因的长度进行校正，之后再用校正后的这个基因read数与校正后的这个样本的所有read数求商。
目前大家认为TPM更准确；
获取到不同的然后将不同样本的表达量合并为我们常见的表达矩阵。

接着就是常规的分析流程：PCA、聚类、基因差异表达分析、功能富集分析。常规套路。

后续使用R语言分析：

(1)PCA

数据降维，绘制出前两个主成分二维图；可以查看不同批次样本是否存在批次效应，如果存在批次效应，使用combat包去除；

(2)聚类分析

K-means聚类、层次聚类、NMF等，目的是为了将样本聚为几个大类别，更好地分析；具体原理百度；R包：K-means使用consensusClusterplus；层次聚类使用hcluster;NMF 使用nmf；

(3)差异表达分析

目的是为了获取不同类别间的差异基因，通过统计学检验得到。Deseq2和edgeR

(4)功能富集分析

目的是看差异表达基因和哪些功能相关，使用超几何检验得到；R包：clusterProfiler，网站：webgestalt
不推荐使用DAVID网站；
有样本表达量，还可以使用GSEA进行富集分析：GSEA主要是输入文件格式比较麻烦，然后看懂GSEA结果；

最基本的分析流程、期间会有很多问题出现；
需要把重点放在后续分析的方法：如何与癌症相结合; 通过工具实现我们研究目的，会使用基本工具，更重要的是知道为什么要这样分析，想要得到某个预想的结果，需要用什么方法，需要长期的积累与实践。

RNA-seq详细教程链接：
转录组入门传送门

代码整理:

fastq-dump

ls *.sra|awk '{print "fastq-dump --split-3 "$0}'|sh

Trimmomatic
没有adapter:

java -jar /WORK/paratera_3/software/Trimmomatic-0.33/trimmomatic-0.33.jar SE -threads 2 -phred33 -trimlog Temp/SRR035243.log（日志文件） SRP001843/SRR035243/SRR035243.fastq（处理文件） TrimData/SRR035243.trim.fq.gz（结果文件） LEADING:20 TRAILING:20 SLIDINGWINDOW:4:15 MINLEN:36

有adapter:

java -jar /WORK/paratera_3/software/Trimmomatic-0.33/trimmomatic-0.33.jar SE -threads 2 -phred33 -trimlog Temp_fa/output.log   input.fastq TrimData_fa/ouput.trim.fq.gz ILLUMINACLIP:adapter/input_adapter.fa:2:30:10 LEADING:20 TRAILING:20 SLIDINGWINDOW:4:15 MINLEN:36

前后20，最少36bp

samtools 格式转换、排序、构建索引

samtools view -S /disk/zhw/RNAseq/hisat2/${i}.sam -b> ${i}.bam 
samtools sort ${i}.bam -o ${i}_sorted.bam
samtools index ${i}_sorted.bam

比对工具

bowtie
得加上-S 才有头文件

bowtie -S /disk/database/human/hg38/Gencode/bowtie_index/bowtie_genome SRR577511.fastq SRR577511.sam

建立转录本的索引，后面用给RSEM用bam文件定量

bowtie-build  /disk/database/human/hg19/RSEM_bowtie_index/human_RSEM_bowtie.transcripts.fa human_RSEM_trans

bowtie2
单端：

bowtie2 -q -p 24 -x /WORK/paratera_3/database/mm10/bowtie2Index/GRCm38 -U /WORK/paratera_3/data/zhw/Chip_Seq/0_TF_Histone_data/SRP001843/SRR035241/SRR035241.fastq -S SRR035241.sam

双端：

bowtie2 -q -p 24 -x /WORK/paratera_3/database/mm10/bowtie2Index/GRCm38 -1 /WORK/paratera_3/data/zhw/Chip_Seq/0_TF_Histone_data/SRP015991/SRR579523/SRR579523_1.fastq -2 /WORK/paratera_3/data/zhw/Chip_Seq/0_TF_Histone_data/SRP015991/SRR579523/SRR579523_2.fastq -S SRR579523.sam

-q fastq文件输入，-p 24个线程，-x ,索引文件前缀 单端：-U 输入文件,双端-1 输入文件1 -2 输入文件2 -S 输出文件

STAR
索引构建:

STAR --runThreadN 20 --runMode genomeGenerate --genomeDir  /disk/database/human/hg38/Gencode/STAR_index_150bp --genomeFastaFiles /disk/database/human/hg38/Gencode/genome.fa --sjdbGTFfile /disk/database/human/hg38/Gencode/gencode.v25.annotation.gtf

比对

STAR --runThreadN 3 --outSAMtype BAM SortedByCoordinate  --genomeDir /disk/database/human/hg19/STAR_index   --readFilesCommand zcat --readFilesIn s5_H7YHTCCXY_L6_1.clean.fq.gz s5_H7YHTCCXY_L6_2.clean.fq.gz --outFileNamePrefix s5

--readFilesCommand 输入文件格式，zcat解压
--outSAMtype 输出文件格式BAM
SortedByCoordinate 排好序

hista2

hisat2 -t -x /disk/zhw/RNAseq/index/hg19/genome -1 /disk/zhw/RNAseq/rawdata/SRR1294495_1.clean.fq -2 /disk/zhw/RNAseq/rawdata/SRR1294495_2.clean.fq -S SRR1294495.sam

** tophat2**

tophat -p 8 -G /data/database/human/hg19/gencode.v27lift37.annotation.gtf -o  /data/data/zhw/align/tophat  /data/database/human/hg19/Bowtie2Index/genome  /data/data/zhw/align/fastq/SRR1294495_1.clean.fq /data/data/zhw/align/fastq/SRR1294495_2.clean.fq

RSEM
官方文档：
http://deweylab.github.io/RSEM/rsem-calculate-expression.html

教程：
https://github.com/bli25ucb/RSEM_tutorial#intro

Reference:

bowtie2:

rsem-prepare-reference --gtf gencode.v27lift37.annotation.gtf (注释文件) --bowtie2 --bowtie2-path  /disk/soft/bowtie2-2.3.2 (bowtie2路径)  /disk/database/human/hg19/RSEM_bowtie2/hg19.fa (参考基因组)  /disk/database/human/hg19/RSEM_bowtie2/ref/human_0 (保存路径)

bowtie：

rsem-prepare-reference --gtf /disk/database/human/hg19/gencode.v27lift37.annotation.gtf --bowtie --bowtie-path  /disk/soft/bowtie-0.12.7 /disk/database/human/hg19/RSEM_bowtie2/hg19.fa /disk/database/human/hg19/RSEM_bowtie_index/human_RSEM_bowtie

STAR的参考要加上

rsem-prepare-reference --gtf /disk/database/human/hg38/Gencode/gencode.v25.annotation.gtf --star --star-path /disk/soft/STAR -p 20 /disk/database/human/hg38/Gencode/genome.fa /disk/database/human/hg38/Gencode/RSEM_STAR_150bp/human_150bp

比对加定量：

bowtie:

rsem-calculate-expression -p 20 --paired-end /disk/zhw/cross_talk/GSE57872/testfastqfile/SRR1294493_1.good.fq  /disk/zhw/cross_talk/GSE57872/testfastqfile/SRR1294493_2.good.fq /disk/database/human/hg19/RSEM_bowtie_index/human_RSEM_bowtie ./SRR1294493

设置botwtie2,索引应该的bowtie2的参考索引文件，因为会调用bowtie2

rsem-calculate-expression -p 20 --paired-end --bowtie2 --bowtie2-path /disk/soft/bowtie2-2.3.2 --output-genome-bam --seed-length 15  /disk/zhw/cross_talk/GSE57872/afterqc/SRR1294809_1.good.fq  /disk/zhw/cross_talk/GSE57872/afterqc/SRR1294809_2.good.fq  /disk/database/human/hg19/RSEM_bowtie2/ref/human_0 ./SRR1294809

index 加上文件前缀

STAR:
--star-gzipped-read-file 输入文件为压缩文件，解压
--star-bzipped-read-file
--star-path /disk/soft/STAR ：STAR所在文件夹名称
/disk/database/human/hg19/RSEM_STAR/GRCh37：GRCh37是文件前缀

rsem-calculate-expression -p 20 --paired-end --star  --star-path /disk/soft/STAR /disk/zhw/cross_talk/GSE57872/testfastqfile/SRR1294493_1.good.fq  /disk/zhw/cross_talk/GSE57872/testfastqfile/SRR1294493_2.good.fq /disk/database/human/hg19/RSEM_STAR/GRCh37 ./SRR1294493

输入bam,直接定量：

rsem-calculate-expression -alignments -p 6 --estimate-rspd /disk/zhw/cross_talk/GSE57872/bowtie/SRR1294492.bam /disk/database/human/hg19/RSEM_bowtie_index/human_RSEM_bowtie

kallisto
https://pachterlab.github.io/kallisto/starting

构建index

kallisto index -i transcripts.idx /disk/soft/kallisto_linux-v0.43.1/test/transcripts.fasta.gz
kallisto index -i hg38_transcripts.idx /disk/database/human/hg38/Ensembl/Homo_sapiens.GRCh38.cdna.all.fa

定量

kallisto quant -i transcripts.idx -o output -b 100 /disk/soft/kallisto_linux-v0.43.1/test/reads_1.fastq.gz /disk/soft/kallisto_linux-v0.43.1/test/reads_2.fastq.gz

七、ChIP-seq流程分析

ChIP-seq原理：具体请百度，大致是将需要结合的蛋白质与序列结合后，使用蛋白质抗原抗体结合技术，将那部分结合了蛋白质的序列拉下来，进行测序；

RNA-seq 前面主要是为了获取表达量；而ChIp-seq分析是为了研究某个位点有蛋白质结合以及峰值；
分析流程前面一样，都是先对原始数据进行质量控制，查看数据质量情况，如果有比较差的质量的read,需要进行修剪，接着就是将read比对到基因组上；
RNA-seq后面则是进行转录本的计数，即定量；而ChIP-seq一般是通过macs2寻找Peaks富集区；

ChIP-seq 详细教程链接：
转录组入门传送门

RNA-seq和ChIP 对比，注意，这里刚好是做super enhancer,且主要是做差异peak；

代码整理:
merge重复样本 **

samtools merge NKTL010_S26.sam NKTL010_S26_L005_R1_001.sam NKTL010_S26_L006_R1_001.sam

callpeak

macs2 callpeak -t NKTL011_S27.sam -c NKTL010_S26.sam -f SAM -g hs -n NKTL011_S27 -B -q 0.01 >NKTL011_S27.macs2.log &

ROSE 寻找超级增强子

cd /disk/soft/rose
python /disk/soft/rose/ROSE_main.py -g hg38 -i /disk/zhw/zuo/chipseq/NKTL_ChIPseq_data/NKTL011_S27_peaks.gff -r /disk/zhw/zuo/chipseq/NKTL_ChIPseq_data/sam/bam/sortedbam/NKTL011_S27_L006_R1_001.sorted.bam -c /disk/zhw/zuo/chipseq/NKTL_ChIPseq_data/sam/bam/sortedbam/NKTL011_S27_L005_R1_001.sorted.bam -o /disk/zhw/zuo/chipseq/NKTL_ChIPseq_data/sam/bam/NKTL011_S27 -s 12500 >NKTL011_S27.log

mergebed 做intersect

mergeBed -i SNK6_NKTL011_Enhancers_sorted.bed -c 4 -o collapse> SNK6_NKTL011_Enhancers_mergeBed.bed

注释到附近基因

closest-features --closest Enhancers_SNK6_NKTL011.bed /disk/database/human/hg38/Gencode/gencode.v25.annotation.sorted.bed>Enhancers_SNK6_NKTL011_closegene.bed

排序，提取，替换

sortBed -i  NKTL002_S18.sorted.bed |awk '{print $1,$2,$3,$6}' |sed 's/ /	/g'>NKTL002_S18.sorted.read.bed

MAnorm寻找差异peak

sh ./MAnorm.sh /disk/zhw/zuo/chipseq/NKTL_ChIPseq_data/mergebed/NKTL008_S24_peaks_Gateway_Enhancers_sorted_peak.bed /disk/zhw/zuo/chipseq/NKTL_ChIPseq_data/mergebed/NKTL002_S18_peaks_Gateway_Enhancers_sorted_peak.bed NKTL008_S24.sorted.read.bed  NKTL002_S18.sorted.read.bed 150 148

bam文件转bed文件

bamToBed -i NKTL011_S27.sorted.bam >NKTL011_S27.sorted.bed

转换bam文件格式为tdf

/disk/soft/IGVTools/igvtools count  NKTL011_S27_L006_R1_001.sorted.bam NKTL011_S27_L006_R1_001.sorted.bam.tdf hg38

八、芯片数据分析

芯片测序原理：
与RNA-seq不同，芯片主要是通过结合在cDNA上的荧光或者同位素的信号强度来代表表达量；

预处理：
主要是背景处理、归一化和汇总，分别都有不同的函数，但是也有三合一的，主要是rma和mas5；

背景处理（background adjustment）: MAS和RMA
归一化处理（normalization，或称为“标准化处理”）
汇总（summarization）

rma函数和ma5函数 三合一:
rma:归一化处理使用分位数法，而汇总方法使用medianpolish
Mas5: expresso函数和mas方法
具体教程：基因芯片（Affymetrix）分析2：芯片数据预处理

不同平台处理方式：
Affymatix的主要是使用affy包，使用ReadAffy读取文件，rma一步预处理获取表达矩阵；Illumina平台的数据主要是使用limma包，先read.ilmn读取文件，运用rma进行背景标准化：（backgroundCorrect），运用quantile进行组间标准化（normalizeBetweenArrays）
Agilent与Illumina相似，主要是读取文件使用的是read.maimages这个函数；背景标准化和组间标准化函数的一样的，但是单双通道的数据处理方式不一样；而由于芯片数据是使用probeid标记不同的基因，所以需要注释到不同的基因名，一般GEO的不同平台的GPL页面都会有对应的注释文件，把id和芯片数据对应起来就好，也可以去官网找；
获取到表达矩阵后，后续可以进行PCA，聚类和差异表达分析等常规套路；差异表达一般使用limma包

九、全外显子测序数据分析pipeline

总体流程：

质控——比对——merge——去重——重比对——BQSR——位点检测——filteration——annova注释

质控

下机数据(fastq)

第一行:“@”开头+ Illumina 测序标识别符和描述文字，即 read ID
第二行:碱基序列
第三行:以“+”开头，随后为 Illumina 测序标识别符,为了节省存储，一般为空
第四行：对应碱基的测序质量=每个字符对应的 ASCII 值-33（phred33,若phred64需要转格式）

如何判断自己的数据质量分数是phred33还是phred64呢
如果有2个以上的质量字符ASCII值小于等于58，同时没有任何质量字符的ASCII值大于等于75，即判断是Phred+33
如果有2个以上的质量字符ASCII值大于等于75，同时没有任何质量字符的ASCII值小于等于58，即判断是Phred+64

数据质控

1.可用fastqc来看原始reads的质量
2.过滤步骤：
①去除adptor的reads
② 去除 N 的比例大于 10%的 reads
③去除低质量 reads（质量值 Q <= 3 的碱基数占整个 read 的 50％以上）
工具：Trimmomatic ,fastx_toolkit
3.测序深度：测序得到的总碱基数于待测基因组大小的比值
覆盖度：测序产出的数据覆盖到目标基因组上的比例

比对

用比对软件将clean reads 与参考基因组做比对生成sam文件 ,并且转为bam格式的文件
工具：** BWA(mem) **

bwa mem -t  -k  -m -R   <refernce> <read1.fa> <read2.fa>  |  samtools view -bS > bam

-t 线程 -k seed中的最大编辑距离 -M 与picard兼容 -R 为reads加头 -b 输出的为 bam文件

bam文件中同一染色体对应的条目按照坐标顺序从小到大进行排序
工具：picard (sortsam)

java -jar picard SortSam I=bam O=sorted.bam SORT_ORDER=coordinate

merge

一个样本可能有多条line,有多个bam文件，因此需要合并为一个bam文件
工具：picard merge
INPUT_FOR_picard_MarkDuplicates+=INPUT="sorted.bam"

去重

在制备文库的过程，会进行PCR反应，但是这些序列仍会比对到参考基因组上，不能作为检测变异的证据，需要通过picard 给那些重复序列加上标签
工具：picard Duplicates Marking

java -jar picard MarkDuplicates REMOVE_DUPLICATES=false MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000  $ INPUT_FOR_picard_MarkDuplicates  OUTPUT=dedup.sorted.bam M=rmdup.metrics.txt

重比对

在indel附近的比对会出现大量的碱基错配，会被误认为是snp,因此对indel附近进行重比对
工具：GATK
①用RealignerTargetCreator根据indel相关数据库来确定进行重比对的区间

java -jar GATK -R reference -T RealignerTargetCreator -known Mills_and_1000G_gold_standard.indels.hg19.sites.vcf -known 1000G_phase1.indels.hg19.sites.vcf -I  dedup.sorted.bam -o realn.dedup.sorted.intervals

②用IndelRealigner在这些区域内进行重新比对

java -jar GATK -R reference -T IndelRealigner -targetIntervals realn.dedup.sorted.intervals -known Mills_and_1000G_gold_standard.indels.hg19.sites.vcf -known 1000G_phase1.indels.hg19.sites.vcf  -I dedup.sorted.bam -o realn.dedup.sorted.bam

BQSR

碱基质校正：bam文件里reads的碱基质量值进行重新校正，使最后输出的bam文件中reads中碱基的质量值能够更加接近真实的与参考基因组之间错配的概率

工具：GATK BaseRecalibrator AnalyzeCovariates PrintReads

① 用BaseRecalibrator根据known sites来生成一个校正质量所需的文件

java -jar GATK  -R reference -T BaseRecalibrator -knownSites dbsnp_137.hg19.vcf   -knownSites  Mills_and_1000G_gold_standard.indels.hg19.vcf  -knownSites 1000G_phase1.indels.hg19.vcf -I realn.dedup.sorted.bam -BQSR before.grp

② 用校正前的文件来生成一个校正后的数据文件

java -jar GATK -R index -T BaseRecalibrator -knownSite 1000G_phase1.indels.hg19.sites.vcf -knownSites Mills_and_1000G_gold_standard.indels.hg19.sites.vcf -knownSite dbsnp_138.hg19.vcf-I realn.dedup.sorted.bam -BQSR before.grp - o post.grp

③ 用AnalyzeCovariates评估前后碱基质量值的比较结果，表格或图形

java -jar GATK -R reference -T AnalyzeCovariates -before before.grp -after post.grp -csv quality.csv -plots quality.pdf

④用PrintReads将经过质量值校正的数据输出到新的bam文件中

java -jar GATK -R refernce -T PrintReads -I realn.dedup.sorted.bam -BQSR before.grp -o checked.bamsampleSub

位点检测

call 变异位点

工具：** GATK **

germline:

java  -jar GATK -T HaplotypeCaller -R inference -I checked.bam --genotyping_mode DISCOVERY -stand_call_conf 30 -o raw.variants.vcf -nct 20  

somatic:

java -jar GATK -T MuTect2 -R inference -I:tumor tumor.bam  -I:normal normal.bam
-pon MuTect2_normal_pon.vcf.gz -o tumor _normal.vcf

filteration

两种过滤不可信位点方法
位点矫正：** VQSR：GATK **推荐样本数目要达到30个以上

去除不可信：Hard filteration：样本数目少，直接控制条件

这里主要是call snv:snp和indel

①先从raw中分别取出snp和indel

java -jar ${GATK} -T SelectVariants -R reference -V sample.raw.variants.vcf -selectType SNP -o sample.snps.vcf 
java -jar GATK -T SelectVariants -R reference -V sample.raw.variants.vcf -selectType INDEL -o sample.indels.vcf

②VQSR

SNP

java -jar GenomeAnalysisTK.jar -T VariantRecalibrator -R reference -input  sample.snps.vcf 
-resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap.vcf 
 -resource:omni,known=false,training=true,truth=true,prior=12.0 omni.vcf     
 -resource:1000G,known=false,training=true,truth=false,prior=10.0 1000G.vcf 
 -resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp.vcf 
-an DP -an QD -an FS  -an SOR  -an MQ -an MQRankSum -an ReadPosRankSum 
-an InbreedingCoeff -mode SNP -tranche 100.0 -tranche 99.9 -tranche 99.0 -tranche 90.0 
-recalFile recalibrate_SNP.recal -tranchesFile recalibrate_SNP.tranches 

java -jar GATK -T ApplyRecalibration -R reference -input sample.indels.vcf -mode SNP  
--ts_filter_level 99.0 -recalFile recalibrate_SNP.recal  -tranchesFile recalibrate_SNP.tranches -o sample.re_snp.vcf 

indel

java -jar GATK -T VariantRecalibrator  -R reference -input   sample.indels.vcf
- resource:mills,known=false,training=true,truth=true,prior=12.0  Mills_and_1000G.vcf 
-resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp.b37.vcf  
-an QD -an DP -an FS -an SOR -an MQRankSum -an ReadPosRankSum -an InbreedingCoeff-mode INDEL -tranche 100.0 -tranche 99.9 -tranche 99.0 -tranche 90.0 --maxGaussians 4 
-recalFile recalibrate_INDEL.recal -tranchesFile recalibrate_INDEL.tranches 

java -jar GenomeAnalysisTK.jar  -T ApplyRecalibration   -R reference -input sample.indels.vcf  recalibrated_snps_raw_indels.vcf  -mode INDEL --ts_filter_level 99.0 
-recalFile recalibrate_INDEL.recal  -tranchesFile recalibrate_INDEL.tranches
 -o sample .re_indel.vcf 

而对于体细胞突变的检测

把Normal的位点合并在一起，从而过滤掉在正常组织中也存在的变异位点

java -jar GATK -T CombineVariants -R reference 
-V normal1.vcf -V normal2.vcf ............-v normaln.vcf 
 -minN 2 --setKey "null" --filteredAreUncalled --filteredrecordsmergetype KEEP_IF_ANY_UNFILTERED -o MuTect2_PON.vcf

③hard filteration

snp

java -jar GATK -T VariantFiltration -R reference  -V sample.snps.vcf  
 --filterExpression "QD < 2.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0"   --filterName "my_snp_filter"  -o sample.filtered_snps.vcf

indel

java -jar GATK -T VariantFiltration -R reference.fa -V sample.indels.vcf 
--filterExpression "QD < 2.0 || FS > 200.0 || ReadPosRankSum < -20.0" 
--filterName "my_indel_filter" -o filtered_indels.vcf 

给条件，硬过滤相当于给符合上面条件的加上标签my_snp_filte/my_indel_filter，不符合PASS标签。

解释：

germline
VQSR：根据已有的真实变异位点:HapMap3、Omni 2.5M、 SNP芯片已知的数据去训练高斯混合模型，然后用来评估自己的变异位点，根据聚类结果赋予所有变异位点相应的VQSLOD值。

Hard filteration: INDEL——QD 、FS、ReadPosRankSum
SNP——QD、FS、 MQ 、MQRanksum 、 ReadPosRankSum
QD：单位深度变异质量值，越大可信度越高
FS：只含有变异的read以及只含有参考序列碱基的read是否存在着明显的正负链特异性，值越小越好。
MQ：所有比对至该位点上的read的比对质量值的均方根
保留： PASS 去掉：带上不合格的标签 grep取出pass

ANNOVAR注释

注释软件用table_annovar.pl根据一直的数据库注释call出来的位点

perl annovar/table_annovar.pl  filted.vcf  annovar/humandb/ -buildver hg19 -out annot.variants.vcf -remove -protocol refGene,cosmic70,avsnp147,
ALL.sites.2015_08_edit,EAS.sites.2015_08_edit,esp6500siv2_all,exac03,ljb26_all -operation g,f,f,f,f,f,f,f -nastring . --vcfinput