一起学生信分析 RNA-Seq上游 篇

参考文章：

https://zhuanlan.zhihu.com/p/345896914

RNA-Seq分析介绍

转录组是指某特定细胞类型产生的所有转录本的集合。转录组研究能够从整体水平研究基因功能以及基因结构，揭示特定生物学过程以及疾病发生过程中的分子机理，已广泛应用于基础研究、临床诊断和药物研发等领域。RNA测序(RNAseq)自诞生起就应用于分子生物学，帮助理解各个层面的基因功能。现在的RNA-seq更常用于分析差异基因表达(DGE, differential gene expression)，而从得到差异基因表达矩阵。RNAseq在过往十年里逐渐成为全转录组水平分析差异基因表达和研究mRNA差异剪接必不可少的工具。

分析流程

1. 质量控制 (Quality Control)

   1. FastQC + Trimmomatic + Cutadapt
   2. Fastp

2. 比对到参考基因组 (Mapping to a reference genome)

   1. STAR
   2. HISAT2
   3. Tophat2

3. 基因表达量定量 (Quantification)

   1. featureCounts 或者 HTSeq -- 从比对结果中计数基因的reads。
   2. Salmon 能够直接从FASTQ文件进行基因或转录本的表达量定量，不需要预先的比对步骤。

4. 差异表达分析 (Differential Expression Analysis)

   1. DESeq2
   2. edgeR
   3. limma

5. 功能富集分析 (Functional Enrichment Analysis)

   1. clusterProfiler 或 GOseq -- 进行基因本体(GO)富集分析。
   2. GSEA - 基因集富集分析。

6. 可视化

快速上手RNA-Seq分析流程 -- 上游

STEP01 质量控制 (Quality Control) -- 使用fastp

# 新建 clean 文件夹存放 fastp 质控之后的数据
mkdir -p clean

# 单样本处理实例，默认选项
fastp -w 20 -i raw/sample1_R1.fq.gz -I raw/sample1_R2.fq.gz \
-o clean/sample1_R1.fastq -O clean/sample1_R2.fastq \
--report_title "${i} fastp report" \
--json clean/${i}_fastp.json --html clean/${i}_fastp.html

# 多样本处理
# fastp质控
for i in `tail -n +2 metadata.txt| cut -f1`
do

fastp -w 20 -i raw/${i}_R1.fq.gz -I raw/${i}_R2.fq.gz \
-o clean/${i}_R1.fastq -O clean/${i}_R2.fastq \
--report_title "${i} fastp report" \
--json clean/${i}_fastp.json --html clean/${i}_fastp.html

done

STEP02 比对到参考基因组 (Mapping to a reference genome) -- 使用STAR

比对软件的使用一般分为两步：

1. 建立索引

# --genomeDir 这里是索引存储文件夹  
# --runMode genomeGenerate 任务为建立索引

STAR \
--runThreadN 50 \
--runMode genomeGenerate \
--genomeDir ./index \
--genomeFastaFiles ../genome/GCA_014117465.1_ASM1411746v1_genomic.fna \
--sjdbGTFfile ../genome/genomic.gtf \
--sjdbOverhang 100

2. 比对到基因组

# 单样本
STAR 
--runThreadN 5 #线程数为5
--genomeDir ./index #索引位置
--readFilesCommand cat #读取文件
--readFilesIn clean/sample1_R1.fastq clean/sample1_R2.fastq
#输入质量过滤后的文件
--outFileNamePrefix ./STAR/sample1_
#输出文件路径与命名方式
--outSAMtype BAM #输出BAM格式
SortedByCoordinate #基于位置对输出文件排序
--outBAMsortingThreadN 5 #输出文件排序使用线程数为5
--quantMode TranscriptomeSAM #同时生成基于转录本的比对文件
GeneCounts #计数


# 多样本
for i in `tail -n +2 metadata.txt | cut -f1`
do

STAR \
--runThreadN 40 \
--genomeDir ./index \
--readFilesCommand cat \
--readFilesIn clean/${i}_R1.fastq clean/${i}_R2.fastq \
--outFileNamePrefix ./STAR/${i}_ \
--outSAMtype BAM \
SortedByCoordinate \
--outBAMsortingThreadN 10 \
--quantMode TranscriptomeSAM \
GeneCounts

done

基因表达量定量 (Quantification) -- 使用featureCounts

# -a 指定基因组注释文件
featureCounts \
-a ../genome/genomic.gtf \
-p \
-T 50 \
-o result/featureCounts_InOutput/matrix.txt \
STAR/*.sortedByCoord.out.bam

至此，就可以得到差异表达矩阵，并进行相应的下游分析。