宝塔服务器面板，一键全能部署及管理，送你10850元礼包，点我领取

纯粹就是记录下如何在Window下少写代码完成RNA-seq,并没有详细讲解转录组分析的原理和每一步背后的含义。想深入了解可以翻生信技能树，生信宝典，组学大讲堂等等公众号。本文里的命令都是用powershell执行的，cmd当然也可以。
本文主要参考CJ写的完结 | TBtools RNA-seq 数据分析系列插件
从下载数据开始到完成初步结果大概花了15个小时，大致为11步。
本文是从0开始RNA-seq,本文前面几步是找别人的转录组测序数据，不需要可以跳过前面几步。

工具准备

TBtools
Sra Toolkit
FastQC
Trimmomatic
Kallisto

1、SRA 数据查询与整理

在NCBI 的SRA 中找了下水稻的测序数据，看到水稻测序共有9w+条个，转录组测序是1.7w+。下载了xml查看这些项目的基本信息。用SRA XML to Table 可以转换为可查看的模式，并附带下载链接。

2、SRA 数据下载

找到感兴趣的测序数据就可以下载了，用啥下都行，迅雷也行，就是多个数据的话不太友好。批量下载，用sra toolkit工具也可以，就是速度不太行(可能我这网不好)。有些项目提供SRA和reads的fastq，有些只有SRA。当然Aspera 是最好的选择，可以参考使用Aspera高速下载SRA/ENA测序数据。TBtools也内置了Aspera接口。

3、SRAtoFastq

如果拿到的是SRA文件就需要转为fastq格式才能进行后续分析。常规做法是用NCBI自己的sra toolkit工具(三个操作系统都有)
sra toolkit下载地址
具体可以使用方法可以参考Windows系统下载SRA数据，使用sratoolkit工具
如果是第一次使用需要在cmd里输入命令 激活

vdb-config –interactive

SRA转Fastq 的命令很简单，不过先要看好NCBI上你下的SRA数据是单端测序还是双端测序，在SRA文件目录打开powershell或cmd。双端测序会生成两个文件，单端是一个文件。

fastq-dump –split-3 .\SRR15254739 #双端测序用fastq-dump .\SRR15254739 #单端测序

当然TBtools上也能实现相同功能，需要下载插件SRA to Fastq 。这个插件是要license 的，有需要可以微信号找下WatchBio，找这个人要。

4、FastQC

拿到reads后就是要看下质量。常规就是用FastQC。这个软件先安装java。安装和使用教程参考
windows 安装 fastqc; fastqc
下载链接
打开bat文件，会弹出一个界面就可以开始质检了，是可以一次拖多个文件的。

每个结果的具体含义都有很多教程了，这里就不重复了。
当然TBtools上也能实现相同功能，需要下载插件FastQC 。当然这个插件也是要license 的。

5、数据过滤

数据过滤的软件有很多，这里用的是Trimmomatic，也是基于java开发的。下载地址 解压后就能用，就是没有GUI，要在cmd或powershell里使用。官网已经给出了使用教程,百度也有很多。
学习使用一款数据质控软件（Trimmomatic）
在Trimmomatic文件目录下打开powershell

###使用方法java -jar .\trimmomatic-0.39.jar###双端测序是两个输入，4个输出，接头一般是phred33,具体还是要问公司java -jar trimmomatic-0.39.jar PE -threads 16 -phred33 -trimlog F:\ascp\SRR15254739_trim.log F:\ascp\SRR15254739_1.fastq F:\ascp\SRR15254739_2.fastq F:\ascp\SRR15254739_1_paired.fq.gz F:\ascp\SRR15254739_1_unpaired.fq.gz F:\ascp\SRR15254739_2_paired.fq.gz F:\ascp\SRR15254739_2_unpaired.fq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:8:true SLIDINGWINDOW:5:20 LEADING:3 TRAILING:3 MINLEN:36

也可以拿生成的文件再质控看看
当然TBtools上还是能实现相同功能，需要下载插件Trimmomatic。当然这个插件也是要license 的。
不知道接头的话可以参考Trimmomatic | 点点点，测序原始数据质控，技能√get的方法

6、Kallisto

这里用的是Kallisto来做从reads到表达矩阵。这里参考的是RNA-seq无比对直接定量(Kallisto – sleuth流程)；Kallisto下载地址
用的转录本是从MSU上下的水稻cdna序列 all.cdna
在Kallisto文件目录下打开powershell

###建立引索.\kallisto index -i E:\IRGSP-1.0_representative\rice.idx E:\IRGSP-1.0_representative\all.cdna#双端测序 i是引索目录 t 是线程数 b 是抽样数.\kallisto quant -i E:\IRGSP-1.0_representative\rice.idx -o F:\ascp\ -t 4 -b 100 F:\ascp\SRR15254739_1_paired.fq.gz F:\ascp\SRR15254739_2_paired.fq.gz #单端测序.\kallisto quant -i index -o output –single -l length -s SD file.fq.gz

如果要FPKM的话可以用count来算，TBtools的RPKM/FPKM or TPM 可以计算。基因长度可以用Genome Length Filter 计算。不同标准化介绍可以看RNA-Seq分析|RPKM, FPKM, TPM, 傻傻分不清楚？

然后手动把不同样本的表达矩阵合并起来
用Kallisto的主要原因就是简单，需要的内存小，大部分笔记本都能完成。
常规的策略可以看RNA-seq分析工具最优组合
但Kallisto无法发掘新的转录本,准确度受到注释精度影响。CJ也开发了Hisat2+StringTie插件可以完成常规比对流程。

7、转录本合并

Kallisto得到的是转录本水平，大多数人感兴趣的是基因水平，所以需要把不同转录本合并一下，这个用TBtools的Tran Value Sum 就行了。对应关系需要注释文件，我这用的是all.gff3。提取可以用excel(选择mRNA，去重复，分列)，也可以用GXF Position Extract。具体参考汇总 | 转录本表达矩阵 到 基因表达矩阵。

总而言之，到这一步，基因表达矩阵就算弄完了。接下去就是差异表达分析了。

8、差异表达分析

计算差异表达基因和注释常用的是R包，edgeR、limma、DESeq2 这三个的R包教程非常多，感兴趣可以看下
edgeR、limma、DESeq2三种差异表达包比较（RNA-seq数据）;
一文解决RNA测序资料的差异分析（limma,deseq,edger包,以及没有生物学重复情况下差异分析）。TBtools上也更新了以DESeq2为基础的差异表达基因分析插件(Differential Gene Expression Analysis-DESeq2 Wapper)。
由于本人比较懒，找到个在线分析的网站高颜值绘图工具ImageGP具体参考这篇文章在线差异基因分析 （支持自动和指定批次校正）
点limma DE analysis 输入表达矩阵和分类数据就可以了（需要注册,建议先把数据保存到自己的用户中心里，不然网页容易卡，卡，卡）
由于是基于limma的，所以数据是要有生物学重复的，我就编了点数据增加重复。这个网站还是很强大的，另外还有很多有用的参数，有时间再研究研究。

把差异基因给下载下来就可以做富集分析了，这里共找到3k+个差异表达基因。不得不说，这图还是画的挺好看的，也提供原文件，如果不满意可以另外画。

9、GO

TBtools也能完成GO注释详情请看：零基础-绝对完成GO/KEGG pathway富集分析-和-绘图
背景集能在PlantGSEA上找到
这里我还是偷懒了，用了个在线网站AGRIGO2 支持MSU,挺快的。
还有其他的GO网站
PANTHER，MSU要转Uniprot ID。
GOEAST,也支持MSU，教程可以看去东方，最好用的在线GO富集分析工具。

10、KEGG

TBtools也能完成KEGG分析，和上面流程应该差不多，这里推荐两个在线KEGG的网站：
KOBAS支持多个物种，最多只能输入3000个基因。
g:Profiler也有R包、python库，同时可能输出GO结果。这里就是用的这个，结果文件也能下载，不满意就自己画。也可以进行ID转换。
水稻的话，KEGG号要转成RAP的。因为KEGG上只有日本晴的注释，分别是RAP的和NCBI。水稻ID转换也是很头疼的事儿，NCBI的命名方式太逆天了还很难找到ID对应的文件，建议还是用RAP的。

11、GSEA

都有了GO和KEGG富集了，为啥会有GSEA基因富集分析呢，主要还是因为大家发现在一条通路里基因的表达量有升高有下降，而超几何分布检验只是筛差异大的基因，这样无法很好解释这条通路或功能到底是被抑制还是增强了。
这个软件操作也很简单和GO差不多，也有很多教程。可以看看GSEA-基因富集分析
TBtools上也有对应功能点点点 | 真香！Simple GO GSEA 富集分析 ~
今天太晚了，以后再说吧。

[待填坑]other

补充个 hisat2 + featureCounts 得到表达矩阵的流程，这个是真比对，不是kallisto那种伪比对，能够比对到新的可变剪切。

最后

我发现最耗时的过程还是 sra转fastq的过程，其他操作都挺快的，就是点点点，写三条命令就齐活了。两个样本差不多占了50G的内存，其实最重要的还是表达矩阵，中间文件都可以删了。转录组技术如此成熟，其实在自己笔记本上也完成也不会太费事，当然是在样本少的情况下，我感觉十几个样本的话这样点点点也能接受，样本太多的话还是在服务器上跑省事儿，TBtools也集成了转录组分析的工具，可以在上面完成所有操作，但CJ的意思貌似不会在主体中集成这些插件，有点可惜。最后的最后，转录组分析最重要的还是实验设计，好的实验设计做出的结果才有意义。