使用salmon和sleuth进行小麦RNA-seq差异表达分析
blast大家一定熟悉吧, 最近有大神爆出一个bug。也即参数“-max_target_seqs”并不是你理解的那样。下图红色箭头处是NCBI上的设置选项。根据NCBI的官方文件，这个选项是返回前N个最匹配的结果。如果我们设定为1，则表示返回数据库中最匹配的那个。但实际上却不是，该选项意味着返回数据库中第一个匹配良好的结果（simply returns the first good hit found in the database, not the best hit as one would assume ）。据说该bug最早在2015年被发现，但是迄今NCBI还没有任何改动，哪怕是对“-max_target_seqs”参数所表示的实际意思进行修改。
image-20180926112536008
目前的解决办法就是事后处理，即blast运行完之后再进行筛选。
好了，我们言归正传。上次我们在“评估salmon和kallisto在小麦RNA-seq定量中的异同”中论证了salmon能够充分区分小麦的同源基因，哪怕只有一个SNP也能够有效区分。同时我们也发现，早期版本的kallisto有重大bug，要尽快升级到最新版本。以前使用kallisto进行定量时，有专门的软件sleuth进行差异表达分析。那么sleuth是何方神圣？sleuth2017年年中被发表在nature methods。谷歌显示至今大约被引用116次。自然要比我们常用的DESeq2，edgeR等软件的准确率要高。
image-20180926115001362
下面我们就就重点谈一谈如何使用这两个软件进行差异基因的表达分析。在实际使用过程中发现，这里边很多坑。特此记录一下。
第一部分 参考转录组
参考转录组也即中国春参考基因组上注释出来的基因序列。参考转录组当然越全面越准确才好。变好总得需要一个过程。目前1.0版本的数据有如下几个不足：首先，基因的3’和5’端比较短，没有到头；第二是，很多基因具有多个转录本，1.0给的比较少；第三，遗漏了不少转录本。
具体到我们手里的数据，我们可以采用以下做法。也即先找出自己数据中的新的转录本和新基因，然后补充进参考转录组进行下面的分析。这样有利于发现，我们数据中独特的转录本或基因。当然了，如果只用中国春1.0的参考转录组，对最终的结果影响也不大，例如，功能富集分析等。这些比参考转录组多出来的基因，这些新基因往往是功能未知基因，没有保守结构域等，可以根据表达等筛选比较有趣的进行具体的研究。
下面，我们谈一谈如何获取相对中国春1.0的新转录本或新基因。
1 使用STAR将reads映射至中国春参考基因组。这一步也可以使用其他mapping软件，如hisat2, bowtie2, bbmap等。因为处理的样本较多，我这里使用python写了一个循环处理。熟悉shell的也可以写shell脚本。
#!/usr/bin/env python# -*- coding: utf-8 -*-__author__ = 'wheatomics'import subprocesswith open('input.txt', 'r') as f: for line in f: line = line.strip().split() fq1, fq2 = line print fq1, fq2 proc = subprocess.Popen( ['STAR', '--twopassMode', 'Basic', '--genomeDir', '/data/genome_whole/', '--runThreadN', '30', '--limitSjdbInsertNsj', '5000000', '--outSAMtype', 'BAM','SortedByCoordinate', '--twopass1readsN', '-1','--readFilesCommand', 'zcat', '--sjdbOverhang', '100', '--outSAMattributes', 'All', '--outSAMmapqUnique', '60', '--outFilterMismatchNmax', '10', '--readFilesIn', fq1, fq2, '--outSAMattrRGline', 'ID:' + fq1.split('/')[-1].split('.')[0], 'SM:' + fq1.split('/')[-1].split('.')[0], 'PL:ILLUMINA','--outFileNamePrefix',fq1.split('/')[-1].split('.')[0]], shell=False) proc.wait()# 注意这个参数 --outSAMmapqUnique 60，这是为了后面使用GATK call variants. 如果没有这个需求，可以不用设置，但是如果call variants又只能将一整条染色体分成两部分。所以不能两全其美。
2 筛选bam文件。这里将bam文件里unmapped的reads去掉；因为，我们的样品来自中国春的叶片，所以大于一个mismatch的reads也去掉；去掉mapping长度小于80的reads;只保留proper pair的reads。
for i in {CS_CT_1_1,CS_CT_2_1,CS_CT_3_1,ABA1h_1_1,ABA1h_2_1,ABA1h_3_1,ABA12h_1_1,ABA12h_2_1,ABA12h_3_1,ABA24h_1_1,ABA24h_2_1,ABA24h_3_1};do bamutils filter $i.sorted.bam $i.filtered.sorted.bam -minlen 80 -mapped -properpair -mismatch 1;done
顺带交代下样本的背景。这里是对中国春叶片外施ABA后，分别在1h,12h,24h取样。每个时间点3个重复。这是我现编的。
3 组装转录本。输入文件是上面筛选过的bam文件。软件使用scallop。scallop是去年发表在nature biotechnology上的一款有参转录本组装软件。
image-20180927145240879
据作者的评测要比过往的软件好。但这并不意味着使用二代测序reads组装出来的转录本就是无比正确的，这一点要铭记，特别是当要针对具体的某条转录本进行研究时。
On 10 human RNA-seq samples, Scallop produces 34.5% and 36.3% more correct multi-exon transcripts than StringTie and TransComb, and respectively identifies 67.5% and 52.3% more lowly expressed transcripts. Scallop achieves higher sensitivity and precision than previous approaches over a wide range of coverage thresholds.
for i in {CS_CT_1_1,CS_CT_2_1,CS_CT_3_1,ABA1h_1_1,ABA1h_2_1,ABA1h_3_1,ABA12h_1_1,ABA12h_2_1,ABA12h_3_1,ABA24h_1_1,ABA24h_2_1,ABA24h_3_1}; do scallop -i ${i}.filtered.sorted.bam -o ${i}.gtf; done
上述命令运行完之后会产生12个gtf文件。下面要做的就是合并这12个gtf文件中的转录本。这里使用 stringtie —merge命令。
4 合并转录本。其中iwgsc_refseqv1.0_HCandLC.gtf来自官方的注释文件。mergelist.txt包含上述12个gtf文件的名字，每行一个，共12行。
stringtie --merge -p 8 -G /data/IWGSCv1.0/iwgsc_refseqv1.0_HCandLC.gtf -o meng_merged.gtf ./mergelist.txt
5 与官方注释信息比较，找出新增的转录本或基因。
#比较gffcompare -o gffall -r /data/IWGSCv1.0/iwgsc_refseqv1.0_HCandLC.gtf all_merged.gtf# gffcompare 命令来自 https://ccb.jhu.edu/software/stringtie/gffcompare.shtml# 分离新转录本gtfcuff puniq gffall.all_merged.gtf.tmap all_merged.gtf /data/IWGSCv1.0_RefSeq_Annotations/iwgsc_refseqv1.0_HCandLC.gtf all_merged_unique.gtf# gtfcuff 命令来自 https://github.com/Kingsford-Group/rnaseqtools
6 获取新增转录本的序列。
gffread all_merged_unique.gtf -g data/IWGSCv1.0/genome.fasta -w all_merged_unique.fasta# gffread 命令来自 http://ccb.jhu.edu/software/stringtie/gff.shtml
7 新转录本和中国春参考转录本合并。
cat data/IWGSCv1.0/iwgsc_refseqv1.0_HCandLC.fasta all_merged_unique.fasta > all_merged.fasta
至此新转录本获取完毕。下面选取几个与参考转录本比较下。
上图显示了与参考基因组不同的转录本
第二部分 基因定量
基因的定量使用salmon这个软件。使用的转录本就是上述我们合并而成的。
# indexsalmon index -p 4 -t all_merged.fasta -i all_merged_salmon_index# 定量for i in {CS_CT_1_1,CS_CT_2_1,CS_CT_3_1,ABA1h_1_1,ABA1h_2_1,ABA1h_3_1,ABA12h_1_1,ABA12h_2_1,ABA12h_3_1,ABA24h_1_1,ABA24h_2_1,ABA24h_3_1};do salmon quant -i Uall_merged_salmon_index -l A -1 ../clean/${i}_1_clean.fq.gz -2 ../clean/${i}_2_clean.fq.gz -p 8 --numBootstraps 100 -o ${i}_quant;done# 参数 --numBootstraps 100 一定要放上# 将生成的*_quant文件夹统一放到salmon_result文件夹下面mkdir salmon_resultmv *_quant salmon_result/
至此，定量部分完成
第三部分 差异表达分析
这里使用sleuth。sleuth本来是为kallisto量身定制的，如果salmon想使用sleuth进行差异表达分析，则需要一个R包（wasabi）进行转换。
1 安装wasabi
# 地址 https://github.com/COMBINE-lab/wasabisource("http://bioconductor.org/biocLite.R")biocLite("devtools") # only if devtools not yet installedbiocLite("COMBINE-lab/wasabi")
或者使用 conda install r-wasabi安装
2 结果转换
library(wasabi)sfdirs <- file.path("salmon_result", c("CS_CT_1_quant","CS_CT_2_quant","CS_CT_3_quant","ABA1h_1_quant","ABA1h_2_quant","ABA1h_3_quant","ABA12h_1_quant","ABA12h_2_quant","ABA12h_3_quant","ABA24h_1_quant","ABA24h_2_quant","ABA24h_3_quant"))prepare_fish_for_sleuth(sfdirs)
3 使用sleuth进行差异表达分析
# 这里需要一个分组信息文件“sample_infor.txt”，如下格式，两列。注意第一列的顺序要按字母顺序排列sample conditionABA1h_1 ABA1hABA1h_2 ABA1hABA1h_3 ABA1hABA12h_1 ABA12hABA12h_2 ABA12hABA12h_3 ABA12hABA24h_1 ABA24hABA24h_2 ABA24hABA24h_3 ABA24hCS_CT_1 CS_CTCS_CT_2 CS_CTCS_CT_3 CS_CT# 为了获取基因水平的差异，需要一个转录本和基因对应的文件“transcript_gene_relation.txt”。文件包括两列，第一列是转录本名字，第二列是基因名字。transcript_id gene_idTraesCS1A01G425600.1 TraesCS1A01G425600TraesCS1A01G425700.1 TraesCS1A01G425700TraesCS1A01G425800.1 TraesCS1A01G425800TraesCS1A01G602400LC.1 TraesCS1A01G602400LCTraesCS1A01G602500LC.1 TraesCS1A01G602500LC
我们这里是按时间点取样，可以按照时序分析，不必进行两两比较。
library(sleuth) #加载包library(splines) #加载包sample_id <- dir(file.path(".", "salmon_result")) #加载salmon_result下的文件夹sample_id #显示有多少样品kal_dirs <- file.path(".", "salmon_result", sample_id)kal_dirss2c <- read.table(file.path(".", "sample_infor.txt"), header = TRUE, stringsAsFactors=FALSE)s2c <- dplyr::select(s2c, sample = sample, condition)s2c <- dplyr::mutate(s2c, path = kal_dirs)s2c #注意查看是否一一对应，这个地方一不小心会出现对应错误。也即为啥要求sample_infor.txt里要按第一列排序的原因。design <- model.matrix(~ -1+factor(c(1,1,1, 2,2,2, 3,3,3, 4,4,4))) #此处的design我不是十分确定是否准确？colnames(design) <- c("CS_CT", "ABA1h","ABA12h","ABA24h") designt2g <- read.table(file.path(".", "transcript_gene_relation.txt"), header = TRUE)t2g <- dplyr::rename(t2g,target_id = transcript_id, gene_id = gene_id)so <- sleuth_prep(s2c, full_model = design, num_cores = 1,target_mapping = t2g,aggregation_column = 'gene_id', extra_bootstrap_summary=TRUE,read_bootstrap_tpm=TRUE, gene_mode = TRUE)# 参数 num_cores = 1 一定要设为1，不然会报错。# 如果读取出现错误，可能是因为R包rhdf5的问题，即wasabi和sleuth依赖的rhdf5不同。我这里采取的策略就是分别安装在不同环境里了。或使用conda install --channel bioconda r-sleuth。将wasabi安装在wasabi环境里，conda create -n wasabi r-wasabi。这个环境与sleuth的环境不同so <- sleuth_fit(so)so <- sleuth_fit(so, formula = ~ 1, fit_name = "reduced")so_lrt <- sleuth_lrt(so, "reduced", "full")models(so_lrt)sleuth_table <- sleuth_results(so_lrt, 'reduced:full', 'lrt', show_all = FALSE)sleuth_significant <- dplyr::filter(sleuth_table, qval <= 0.05)head(sleuth_significant)write.csv(sleuth_table,"CS_ABA_sleuth_gene_level.csv",row.names=TRUE,quote=TRUE)write.csv(sleuth_significant,"CS_ABA_sleuth_significant_gene_level.csv",row.names=TRUE,quote=TRUE)sleuth_matrix <- sleuth_to_matrix(so_lrt, 'obs_norm', 'tpm')head(sleuth_matrix)write.csv(sleuth_matrix,"CS_ABA_sleuth_tpm_norm_gene_level.csv",row.names=TRUE,quote=TRUE)
我们也可以进行两两比较
# 首先要明确的是，是那些样本之间进行两两比较。此处是：CS_CT vs ABA1h, CS_CT vs ABA12h, CS_CT vs ABA24h, ABA1h vs ABA12h,ABA1h vs ABA24h, ABA12h vs ABA24h。共六组library(sleuth) #加载包library(stringr) #加载包sample_id <- dir(file.path(".", "salmon_result")) #加载salmon_result下的文件夹kal_dirs <- file.path(".", "salmon_result", sample_id)s2c <- read.table(file.path(".", "sample_infor.txt"), header = TRUE, stringsAsFactors=FALSE)s2c <- dplyr::select(s2c, sample = sample, condition)s2c <- dplyr::mutate(s2c, path = kal_dirs)s2c #注意查看是否一一对应，这个地方一不小心会出现对应错误。也即为啥要求sample_infor.txt里要按第一列排序的原因。t2g <- read.table(file.path(".", "transcript_gene_relation.txt"), header = TRUE)t2g <- dplyr::rename(t2g,target_id = transcript_id, gene_id = gene_id)#下面的语句进入两两比较a <- list(c("CS_CT","ABA1h"), c("CS_CT","ABA12h"), c("CS_CT","ABA24h"), c("ABA1h","ABA12h"),c("ABA1h","ABA24h"), c("ABA12h","ABA24h"))for (xxx in a) {s2b <- dplyr::filter(s2c, condition == xxx[1] | condition == xxx[2])so <- sleuth_prep(s2b, ~ condition, num_cores = 1,target_mapping = t2g,aggregation_column = 'gene_id', extra_bootstrap_summary=TRUE,read_bootstrap_tpm=TRUE,gene_mode = TRUE)so <- sleuth_fit(so)so <- sleuth_fit(so, formula = ~ 1, fit_name = "reduced")so_lrt <- sleuth_lrt(so, "reduced", "full")models(so_lrt)sleuth_table <- sleuth_results(so_lrt, 'reduced:full', 'lrt', show_all = FALSE)sleuth_significant <- dplyr::filter(sleuth_table, qval <= 0.05)head(sleuth_significant)write.csv(sleuth_table,str_c(xxx[1], "vs", xxx[2],"_sleuth_gene_level.csv"),row.names=TRUE,quote=TRUE)write.csv(sleuth_significant,str_c(xxx[1], "vs", xxx[2],"_sleuth_significant_gene_level.csv"),row.names=TRUE,quote=TRUE)sleuth_matrix <- sleuth_to_matrix(so_lrt, 'obs_norm', 'tpm')head(sleuth_matrix)write.csv(sleuth_matrix, str_c(xxx[1],"vs", xxx[2],"_sleuth_tpm_norm_gene_level.csv"),row.names=TRUE,quote=TRUE) }
sleuth还提供了一些画图函数，有兴趣的可以自行尝试。
PCA
sleuth结果里只提供了q value 来筛选差异表达基因，没有倍数变化的结果。这里可以使用excel自行进行计算。
计算 Fold change（倍数变化），十分简单粗暴的方法，计算方法如下：
E = mean (group1) B=mean (group2)
FC= (E-B) / min (E, B)
基因 A 和基因 B 的平均值之差与两者中较小的比值。选择 2-3 倍的基因作为差异表达基因即可。可以根据具体数据，斟酌使用。
提供一个小麦里的文献作为参考吧。今年Cristobal Uauy在BMC Plant Biology上以通讯作者发表的文章，题目是“Ubiquitin-related genes are differentially expressed in isogenic lines contrasting for pericarp cell size and grain weight in hexaploid wheat”。
BMC Plant Biology
q value < 0.05
小麦研究联盟logo
