RNA-Seq分析过程：从FASTQ到图
2019-03-07 本文已影响0人
白夜无寐
Ref: Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis;
写在开头
很多教程会建议大家使用anaconda或者miniconda来安装某些很难很难安装的玩意，比如htseq这货，还有可能的bcftools这货，拿别的怎么都装不上，有时候下源码重编译也没戏。
但是请千万千万千万注意，不要在anaconda的默认环境（root环境，或者base环境）中安装你生物信息学需要的东西。anaconda的环境列表可以使用conda list来列出，也可以在anaconda nevigator的environmen列表中查看，但是无论哪一种，为了你的系统安全着想，请千万在开始任何工作之前，新建一个虚拟环境，这样可以将做生物信息的各种乱七八糟的软件、包和系统隔离开，以避免污染系统PATH的情况，出现这种情况会非常的麻烦，尤其是对os x或者linux系统不熟悉的新手，一旦不新建虚拟环境，path一旦污染，会导致原本已经添加到系统path中的链接无法使用，而同时，你用conda安装的软件也用不了。
另外，流程中的软件/数据路径请务必根据自己的实际情况更改！！
流程千万条，环境第一条；
所以，在开始工作之前，请务必在anaconda中新建一个虚拟环境！！！
利用anaconda安装的时候会遇到很多问题，就个人而言遇到的最大的问题就是samtools......
这个不起眼的工具浪费掉了我一整个下午，最后还是在技能树老师的帮助下才完成了。
最终原因是动态库软连接的问题。就这个问题叙述如下：
1.在使用anaconda之前，我的电脑中已经以源码编译的方式安装了一个samtools。
2.anaconda新建了虚拟环境，然后在虚拟环境中利用conda也安装了一个samtools。
3.如果二者不是同时安装的，那么其实都是可以使用的，在terminal下源码编译安装的话在全局都能用，而单独在conda中安装的话在启动了虚拟环境之后也是没有问题的。
4.如果二者同时安装，无论谁先谁后，问题就来了......
5.解决方案：
如果是在mac系统的话，比较好解决，利用homebrew就可以了。
brew update
brew install openssl
ln -s /usr/local/opt/openssl/lib/libcrypto.1.0.0.dylib /usr/local/lib/
ln -s /usr/local/opt/openssl/lib/libssl.1.0.0.dylib /usr/local/lib/
6.如果是其他发行版的linux的话，也只需要重新安装openssl库就可以了，注意，这里是openssl库而不是openssl的软件，如果是windows的话......抱歉，自己研究去吧......
1.建立索引
1.数据：人类基因组数据（hg19），数据来源于ucsc
2.软件：HISAT2
3.命令：
mkdir hg19
cd hg19
wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz
tar -zxvf chromFa.tar.gz
cat *.fa > hg19.fa
rm chr*.fa
hisat2-build -p 4 hg19.fa genome
#建立基因组索引大约需要半小时左右；平台：MacBook pro 2017；15寸顶配版；cpu为4核8线程；
#建立后会在hg19的文件夹内输出一系列后缀为”.ht2"的文本文件，这是HISAT2的索引文件；
#建立基因组索引是为了加快比对速度；
#然而HISAT2还支持已通过构建snp索引和外显子索引（也即转录组索引），这需要基因组注释文件，然后使用下命令：
extract_exons.py Homo_sapiens.GRCh37.87.gtf > genome.exon
extract_splice_sites.py Homo_sapiens.GRCh37.87.gtf > genome.ss
hisat2_extract_snps_haplotypes_UCSC.py snp151Common.txt > genome.snp
#SNP文件同样来自于ucsc，注意应该是dbSNP文件，这一点来自于HISAT2的官网说明；python脚本已经自带。
#最后统一建成基因组，snp，转录组都有的索引；
hisat2-build -p 4 hg19.fa --snp genome.snp --ss genome.ss --exon genome.exon genome_snp_trans
2.序列比对
1.数据：RNA-Seq测序数据，一般有.sra格式和.fastq格式两种，事实上两种格式区别并不那么大，.sra格式是NCBI由.fastq衍生而来，多了一点描述文本而已，而且用起来还得用sratools转成fastq文件才能使用，所以干脆直接下载fastq用算了。
2软件：比对用的软件很对多，包括TopHat,bowtie2,bwa等，这里由HISAT2一同完成;目前HISAT2是被普遍建议使用的软件，因为速度比较快；
3.命令：
./hisat2 -t -x /Users/coldzero/Desktop/ncbi/hg19/genome_trans -1 /Users/coldzero/Desktop/ncbi/data/human/SRR3589957_1.fastq -2 /Users/coldzero/Desktop/ncbi/data/human/SRR3589957_2.fastq -S /Users/coldzero/Desktop/ncbi/data/human/outdata/SRR3589957.sam
#有两份fastq文件的原因是测序是双端测序，会产生两条fastq文件，当然如果是单端测序那就只有一条了,这里有几个参数的含义为,-t：记录并显示时间；-x：参考基因组；-1:第一条序列文件；-2:第二条序列文件；—S:以sam格式准备输出；其实还有一个隐藏的默认参数-p，紧跟在-t之后，它表示输入的数据是fastq格式，这是HISAT2的默认设置，如果想使用fasta的格式，需要手动将参数改为-f;
#当然了，请注意，序列比对是一个极其消耗内存和cup的步骤，如果电脑不好请千万不要尝试......不然陷入假死等了一夜没结果本人不对此负责。
---------------------------------------------------------------------------------------------
#这是主要的输出结果，以下结果请不要复制
Time loading forward index: 00:00:02
Time loading reference: 00:00:01
Multiseed full-index search: 00:09:08
29720636 reads; of these:
29720636 (100.00%) were paired; of these:
1919822 (6.46%) aligned concordantly 0 times
25504689 (85.81%) aligned concordantly exactly 1 time
2296125 (7.73%) aligned concordantly >1 times
----
1919822 pairs aligned concordantly 0 times; of these:
61657 (3.21%) aligned discordantly 1 time
----
1858165 pairs aligned 0 times concordantly or discordantly; of these:
3716330 mates make up the pairs; of these:
2292119 (61.68%) aligned 0 times
1195950 (32.18%) aligned exactly 1 time
228261 (6.14%) aligned >1 times
96.14% overall alignment rate
Time searching: 00:09:10
Overall time: 00:09:12
3.sam文件操作
1.数据：上一步比对得到的sam结果文件；注意千万要检查上一步的sam文件确实是生成了了。
2.软件：samtools
3.命令：
samtools view -bS SRR3589956.sam >SRR3589956.bam
#这一步将sam文件转换为bam文件，bam文件的好处在于：1.运行速度快，bam文件是二进制文件；2.bam文件进行排序后便于寻找变异；
samtools sort -n SRR3589956.bam -o SRR3589956.sorted.bam
#将排序后的bam文件用于生成bcf文件,排序后的bam文件可用于下游各种分析。
4.转录本组装
1.转录本组装使用的是stringTie软件，该软件可以免费获得；
2.stringTie的输入需要的是刚才sort过的bam文件，需要bam文件是因为这样可以加快运行速度；
3.stringTie的基本参数如下：
4.这一步未必需要进行，根据实际进行取舍。
stringtie <alignen_sorted.bam> [options] *
-h/--help 帮助信息
-v 打开详细模式，打印程序处理的详细信息。
-o [<path/>]<out.gtf> 设置StringTie组装转录本的输出GTF文件的路径和文件名。此处可指定完整路径，在这种情况下，将根据需要创建目录。默认情况下，StringTie将GTF写入标准输出。
-p <int> 指定组装转录本的线程数（CPU）。默认值是1
-G <ref_ann.gff> 使用参考注释基因文件指导组装过程，格式GTF/GFF3。输出文件中既包含已知表达的转录本，也包含新的转录本。选项-B，-b，-e，-C需要此选项（详情如下）
--rf 链特异性建库方式：fr-firststrand(最常用的是dUTP测序方式，其他有NSR，NNSR).
--fr 链特异性建库方式：fr-secondstrand(如 Ligation,Standard SOLiD).
-l <label> 将<label>设置为输出转录本名称的前缀。默认：STRG
-f <0.0-1.0> 将预测转录本的最低isoform的丰度设定为在给定基因座处组装的丰度最高的转录本的一部分。较低丰度的转录物通常是经加工的转录本的不完全剪接前体的artifacts。默认值为0.1。
-m <int> 设置预测的转录本所允许的最小长度.默认值为200
-A <gene_abund.tab> 输出基因丰度的文件（制表符分隔格式）
-C <cov_refs.gtf> 输出所有转录本对应的reads覆盖度的文件，此处的转录本是指参考注释基因文件中提供的转录本。(需要参数 -G).
-a <int> Junctions that do not have spliced reads that align across them with at least this amount of bases on both sides are filtered out. Default: 10
-j <float> 连接点的覆盖度，即设置至少有这么多的spliced reads 比对到连接点(align across a junction)。 这个数字可以是分数, 因为有些reads可以比对到多个地方。 当一个read 比对到 n 个地方是，则此处连接点的覆盖度为1/n 。默认值为1。
-t 该参数禁止修剪组装的转录本的末端。默认情况下，StringTie会根据组装的转录本的覆盖率的突然下降来调整预测的转录本的开始和/或停止坐标。
-c <float> 设置预测转录本所允许的最小read 覆盖度。 当一个转录本的覆盖度低于阈值，则输出文件中不含该转录本。默认值为 2.5
-g <int> 设置ga最小值。 Reads that are mapped closer than this distance are merged together in the same processing bundle. Default: 50 (bp)
-B 应用该选项，则会输出Ballgown输入表文件（* .ctab），其中包含用-G选项给出的参考转录本的覆盖率数据。（有关这些文件的说明，请参阅Ballgown文档。）
如果选项-o 给出输出转录文件的完整路径，则* .ctab文件与输出GTF文件在相同的目录下。
-b <path> 指定 *.ctab 文件的输出路径, 而非由-o选项指定的目录。
注意: 建议在使用-B/-b选项中同时使用-e选项，除非StringTie GTF输出文件中仍需要新的转录本。
-e 限制reads比对的处理，仅估计和输出与用-G选项给出的参考转录本匹配的组装转录本。使用该选项，则会跳过处理与参考转录本不匹配的组装转录本，这将大大的提升了处理速度。
-M <0.0-1.0> 设定。默认值为0.95.
-x <seqid_list> 忽略所有比对到指定的参考序列上的reads，因此这部分的reads不需要组装转录本。 参数 <seqid_list>可以是单个参考序列名称 (如： -x chrM)，也可以是逗号分隔的序列名称列表 (如： -x 'chrM,chrX,chrY')。这可以加快StringTie的组装分析的速度，特别是在排除线粒体基因组的情况下，在某些情况下，线粒体的基因可能具有非常高的覆盖率，但是它们对于特定的RNA-Seq分析可能不感兴趣的。
--merge 转录本合并模式。 在合并模式下，StringTie将所有样品的GTF/GFF文件列表作为输入，并将这些转录本合并/组装成非冗余的转录本集合。这种模式被用于新的差异分析流程中，用以生成一个跨多个RNA-Seq样品的全局的、统一的转录本。
如果提供了-G选项（参考注释基因组文件），则StringTie将从输入的GTF文件中将参考转录本组装到transfrags中。(个人理解：transfrags可能指的是拼接成更大的转录本片段，tanscript fragments)
在此模式下可以使用以下附加选项：
-G <guide_gff> 参考注释基因组文件(GTF/GFF3)
-o <out_gtf> 指定输出合并的GTF文件的路径和名称 (默认值：标准输出)
-m <min_len> 合并文件中，指定允许最小输入转录本的长度 (默认值: 50)
-c <min_cov> 合并文件中，指定允许最低输入转录本的覆盖度(默认值: 0)
-F <min_fpkm> 合并文件中，指定允许最低输入转录本的FPKM值 (默认值: 0)
-T <min_tpm> 合并文件中，指定允许最低输入转录本的TPM值 (默认值: 0)
-f <min_iso> minimum isoform fraction (默认值: 0.01)
-i 合并后，保留含retained introns的转录本 (默认值: 除非有强有力的证据，否则不予保留)
-l <label> 输出转录本的名称前缀 (默认值: MSTRG)
5统计表达量
1.转录本这里用的是featureCounts这是属于subread套件里面的一个组件，该套件中还有subjunc这个神器。
2.同样的工作HTseq也能胜任，但是这货实在太慢了，而且，根据结果来看，二者没有什么区别，而featureCounts这个软件实在是太快了，基本上统计人的只需要30s左右。
3.统计结果就是我们最常见的表达矩阵啦，但这是单样本的表达矩阵，如果要进行后续分析，还要进行转录本信息的合并。
4.转录本技计数并不同于表达量，需要进行标准化才能够真实地展现表达量，
5.命令：
featureCounts -T 4 -p -t exon -g gene_id -a /Users/coldzero/Desktop/ncbi/hg19/Homo_sapiens.GRCh37.87.gtf -o SRR3589957_featurecount.txt SRR3589957.sorted.bam
---------------------------------------------------------------------------------------------
#结果信息如下：
========== _____ _ _ ____ _____ ______ _____
===== / ____| | | | _ \| __ \| ____| /\ | __ \
===== | (___ | | | | |_) | |__) | |__ / \ | | | |
==== \___ \| | | | _ <| _ /| __| / /\ \ | | | |
==== ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
========== |_____/ \____/|____/|_| \_\______/_/ \_\_____/
v1.6.3
//========================== featureCounts setting ===========================\\
|| ||
|| Input files : 1 BAM file ||
|| P SRR3589957.sorted.bam ||
|| ||
|| Output file : SRR3589957_featurecount.txt ||
|| Summary : SRR3589957_featurecount.txt.summary ||
|| Annotation : Homo_sapiens.GRCh37.87.gtf (GTF) ||
|| Dir for temp files : ./ ||
|| ||
|| Threads : 4 ||
|| Level : meta-feature level ||
|| Paired-end : yes ||
|| Multimapping reads : not counted ||
|| Multi-overlapping reads : not counted ||
|| Min overlapping bases : 1 ||
|| ||
|| Chimeric reads : counted ||
|| Both ends mapped : not required ||
|| ||
\\===================== http://subread.sourceforge.net/ ======================//
//================================= Running ==================================\\
|| ||
|| Load annotation file Homo_sapiens.GRCh37.87.gtf ... ||
|| Features : 1195764 ||
|| Meta-features : 57905 ||
|| Chromosomes/contigs : 61 ||
|| ||
|| Process BAM file SRR3589957.sorted.bam... ||
|| Paired-end reads are included. ||
|| Assign alignments (paired-end) to features... ||
|| ||
|| WARNING: reads from the same pair were found not adjacent to each ||
|| other in the input (due to read sorting by location or ||
|| reporting of multi-mapping read pairs). ||
|| ||
|| Pairing up the read pairs. ||
|| ||
|| Total alignments : 34787174 ||
|| Successfully assigned alignments : 23791900 (68.4%) ||
|| Running time : 0.47 minutes ||
|| ||
|| ||
|| Summary of counting results can be found in file "SRR3589957_featurecount ||
|| .txt.summary" ||
|| ||
\\===================== http://subread.sourceforge.net/ ======================//
---------------------------------------------------------------------------------------------
cut -f 1,7 SRR3589957_featureCounts.txt |grep -v '^#' >SRR3589957.txt
---------------------------------------------------------------------------------------------
#从featureCounts的report中可以看到一些东西，比如这里给出的匹配率等信息。在使用了cut之后，我们已经拿到第一列
#为gene_id，第二列为表达量的文本文件了，这就是单样本的表达矩阵了，依次对所有文件进行批处理，然后我们就能得到各
#自的表达谱了。然后随便你用什么办法将它们合并成一个名为matrix的表达矩阵就行，用excel也好，perl也好，python
#也好，R也好，都可以，这里推荐用R，因为这样可以直接接上后续的步骤，而不需要不断切换。
#为了得到表达量，我们应该使用FPKM来表示这里的绝对表达量，首先因为这是PE的数据，不使用RPKM，然后，因为我们拿不
#到所有的转录本数据，所以也不能用TPM。
#另外，这里需要回原本的实验设计中去看一下样本的设计，从EBI的数据网站里可以看到，这三个293的测序数据实际上分为
#两组，其中SRR3589956是control，而SR3589957和SR3589958是 AKAP95的KD，也即treatment，因此在后续分析
#中，应该先做好样本矩阵。当然，这提醒我们，做RNA-Seq之前一定要注意实验的设计；
---------------------------------------------------------------------------------------------
#下面就是接R中的分析步骤了。
#安装DESeq2包
BiocManager::install("DESeq2", version = "3.8")
setwd('/Users/coldzero/Desktop/ncbi/data/human/outdata/')
library(DESeq2)
#配置全局设置，这一步最好加上，避免出现奇奇怪怪的错误；
options(stringAsFactors = FALSE)
#读取各自的count文件，这里的count文件其实就是之前用featureCounts生成的文件；
control1<-read.table('SRR3589956.txt',sep = '\t',col.names = c('gene_id','control1'))
treat1<-read.table('SRR3589957.txt',sep = '\t',col.names = c('gene_id','treat1'))
treat2<-read.table('SRR3589958.txt',sep = '\t',col.names = c('gene_id','treat2'))
#这里可以用head(control)看一下确认一下数据没有问题。
head(control)
control1<-control1[-1,]
treat1<-treat1[-1,]
treat2<-treat2[-1,]
#说明：查看后果然发现有问题，因为featureCount的结果自己带了head，那么有两种解决方法，一种就是在这里的后面删
#掉自带的head，要么就是在读入的时候不设置head。
#然后融合为一个原始count文件,根据gene_id来merge即可，其中可以嵌套
raw_count<-merge(control1, merge(treat1, treat2, by="gene_id"),by = "gene_id")
head(raw_count)
#这里一定要检查raw_count的情况，在gene_id中建议不要出现奇怪的ensembl编号，比如带小数的ensembl编号，如果有，需要进行去除，使用的方法是将小数直接作为空值，相当于直接去掉小数。
a<- gsub("\\_\\d*", "", raw_count_filt$gene_id)
#将.后面的数字替换为空赋值给
ENSEMBL>ENSEMBL <- gsub("\\.\\d*", "", a)
# 将ENSEMBL重新添加到raw_count_filt1矩阵
row.names(raw_count_filt) <- ENSEMBL
raw_count_filt <- cbind(ENSEMBL,raw_count_filt)
colnames(raw_count_filt)[1] <- c("gene_id")
#一定再次检查！！！！
head(raw_count_filt )
#转换gene_id,ensembl是数据库ensembl数据库的储存名称，不好用，我们常用的是symbol。进行基因ID转换的方法非常非常多，这里不建议使用网页工具，尤其是DAVID，因为它的数据库已经茫茫多年没更新了......目前这样的网站简直多如牛毛，这里不再介绍，我们这里用ensembl的接口biomaRt,这是由ensembl官方编写的R的package，更新速度有保证而且结果比较可靠.......或者也可以用org.hs.db来注释，或者clusterProfile注释都可以，都是泛用而且可靠的工具。
library(biomaRt)
library(curl)
#调用boomert接口并使用人类的数据集将ensembl_id转换为symbol。
mart<-useDataset('hsapiens_gene_ensembl',useMart('ensembl'))
#这样会得到一个mart对象，里面储存了人的ensembl的gene信息；
#下面会出现一个很奇怪的问题，这里如果用下面注释掉的命令，得到的my_ensembl_id的变量是错误的，不确定这个问题到底是怎么引起的，按理来说应该不会有这种情况，因此我们换一种写法。
#my_ensembl_id<-row.names(raw_count)（不要用这一个，或者你可以自己尝试一下看看有没有问题）
my_ensembl_id<-raw_count[,1]
options(timeout = 4000000)
hg_symbols<- getBM(attributes=c('ensembl_gene_id','hgnc_symbol',"chromosome_name", "start_position","end_position", "band"), filters= 'ensembl_gene_id', values = my_ensembl_id, mart = mart)
readcount <- merge(raw_count, hg_symbols, by = 'ensembl_gene_id')
head(readcount)
#输出count矩阵
write.csv(readcount,file = 'readcount_all.csv')
#准备样本矩阵这一步请务必自己尝试，一定要弄懂怎么设置。简而言之就是，需要一个单纯的表达矩阵，不带gene_id这一列的，然后将对应的样本的属性（control或者treat）因子化，构成样本信息矩阵，这里需要根据实际的实验设计来调整，也是开始差异分析前的最后一步，务必做对，中间该检查的一定要检查。
#这里一定要弄清楚需要哪些列不需要哪些列，一定要弄清楚！！
row.names(raw_count)<-raw_count$ensembl_gene_id
head(raw_count)
readcount<-raw_count[,-1]
write.csv(readcount,file = 'readcount.csv')
head(readcount)
head(raw_count)
#构建样本信息矩阵，这一步请务必自己尝试，一定要弄懂怎么设置。
condition<-factor(c('control','treat','treat'))
condition
colData<-data.frame(row.names = colnames(readcount),condition)
#开始差异分析，根据DESeq的文档，应当按照标准流程来，先构建dds对象。
#构建dds的时候就会发现，下面的这个函数可以直接食用HTseq的上有数据，当然我们这里是用的自行构建的矩阵，因此选择“frommatrix"
dds<-DESeqDataSetFromMatrix(countData = readcount,colData = colData,design = ~condition)
#这里遇到了一个很迷的错误，提示我作出的mycount每一列的数据类型居然不是数值型，回去检查了一下数据，果然发现，从最开始的文件，read的数据类型就是
#character，这里个人觉得有可能是因为我是用excel处理了一下最开始的txt文件的原因，重新设置一下果然好了......真的是到处都是坑啊。
nrow(dds)
#对表达量极低的基因过滤一下，这里选取的标准是所有样品的基因计数小于1，这也就意味着这个基因在三个样品中read之和都没有1，其实也就相当于都是0。
#在bulk测序中可以将其过滤掉，但是在scRNA-Seq中0值怎么处理就很麻烦了，不能简单�地过滤掉。
dds_filter<-dds[rowSums(counts(dds))>1,]
nrow(dds_filter)
#利用DESeq进行差异分析，分析流程为：
#estimating size factors
#estimating dispersions
#gene-wise dispersion estimates
#mean-dispersion relationship
#final dispersion estimates
#fitting model and testing
dds_out<-DESeq(dds_filter)
#查看分析结果：
res<-results(dds_out)
summary(res)
#结果为：
#out of 25733 with nonzero total read count
#adjusted p-value < 0.1
#LFC > 0 (up) : 1088, 4.2%
#LFC < 0 (down) : 1080, 4.2%
#outliers [1] : 0, 0%
#low counts [2] : 12473, 48%
#(mean count < 40)
#提取分析的结果，我希望得到的是foldchange>2,adjp<0.05（或者0.01）的结果。
table(res$padj<0.05)
res_deseq<-res[order(res$padj),]
summary(res_deseq)
diff_gene_deseq2 <- subset(res_deseq, padj<0.05 & (log2FoldChange > 2 | log2FoldChange < -2))
dim(diff_gene_deseq2)
summary(diff_gene_deseq2)
#提取符合标准的差异表达基因的名称
diff_gene_deseq2 <- row.names(diff_gene_deseq2)
res_diff_data <- merge(as.data.frame(res),as.data.frame(counts(dds_out,normalize=TRUE)),by="row.names",sort=FALSE)
#得到差异分析表格
write.csv(res_diff_data,file = 'res_diff_data.csv',row.names = F)
library(genefilter)
library(pheatmap)
rld <- rlogTransformation(dds_out,blind = F)
write.csv(assay(rld),file="res_diff_data_heatmap.csv")
topVarGene <- head(order(rowVars(assay(rld)),decreasing = TRUE),30)
mat <- assay(rld)[ topVarGene, ]
#减去一个平均值，让数值更加集中。
mat <- mat - rowMeans(mat)
anno <- as.data.frame(colData(rld)[,c("condition","sizeFactor")])
pheatmap(mat, annotation_col = anno)
#至此，我们就拿到了整个流程的第一个初步结果，一张heatmap，从结果来看还可以，聚类也没有什么问题，基因表达矩阵也已经保存，后续高级分析会在下面补充。
最后在此感谢生信菜鸟团和生信技能树团队，主要流程的学习基本都来源于他们，特别感谢技能树的健明老师，帮助解决了途中遇到了无数的坑.......
