本文圍繞RNA-seq學(xué)習(xí)路線進(jìn)行生信入門，主要內(nèi)容有：
☆ RNA-seq方法原理
☆ RNA-seq的生物信息分析
1.數(shù)據(jù)獲取
測序數(shù)據(jù)下載與處理（SRA Toolkit）
測序數(shù)據(jù)質(zhì)控與過濾（fastp）
2.序列比對(duì)（SAMtools、HISAT2）
3.序列組裝（StringTie、TACO）
4.表達(dá)定量和差異表達(dá)分析（Salmon、DESeq2）
5.GO和KEGG富集分析（clusterProfiler）
☆ RNA-seq方法原理
目的是要給mRNA測序，得到樣本的基因表達(dá)信息。
帶Poly（T）探針的磁珠與總RNA進(jìn)行雜交，吸附其中的帶Poly（A）尾巴的mRNA
Mg”離子溶液處理RNA，把RNA打成短片段 被打斷的mRNA片段，用隨機(jī)引物逆轉(zhuǎn)出第一鏈的cDNA，再合成雙鏈cDNA
在雙鏈CDNA的兩端加“A"堿基，并連上"Y“型的接頭
經(jīng)過PCR擴(kuò)增，成為可以上機(jī)的文庫
起始總RNA質(zhì)量控制：用電泳方法。rRNA占有總RNA的大部分，形成的峰越高/尖，RIN（RNA完整度評(píng)分值）越高，8以上質(zhì)量比較好。
測到的RNA片段 mapping到基因組上，進(jìn)行樣品的reads在參考基因上的分布均勻性(Gene coverage)統(tǒng)計(jì)。兩端平衡的時(shí)候表示mRNA降解少（3’高降解多）。
☆ RNA-seq的生物信息分析
一、深度測序數(shù)據(jù)獲取
和EBI、DDBJ組成INSDC，數(shù)據(jù)內(nèi)容相同所以找NCBI就行。
（一）NCBI常用數(shù)據(jù)庫
GenBank：遺傳序列數(shù)據(jù)庫，收集了所有公開的DNA序列及其注釋 GEO (Gene Expression Omnibus)
：收集整理各種表達(dá)芯片數(shù)據(jù)，后來加入了甲基化、lncRNA、miRNA、CNV等其他芯片，還有高通量測序數(shù)據(jù)。文獻(xiàn)中常見GSM和GSE開頭的編號(hào)，分別是GEO
Sample和GEO Series的數(shù)據(jù) PubMed / PMC (PubMed
Central)：前者把測序數(shù)據(jù)和文章聯(lián)系起來，后者可以進(jìn)行全文檢索，無法訪問校園網(wǎng)時(shí)可以替代Web of Knowledge
RefSeq：為所有常見生物提供非冗余、人工挑選過的參考序列，通常包含：參考基因組、參考轉(zhuǎn)錄組、參考蛋白序列、參考SNP信息、參考CNV信息等等
（二）測序數(shù)據(jù)的下載和處理：SRA Toolkit
- 測序數(shù)據(jù)序列格式
（1）FASTA：表示生物序列的文本格式，基因組和EST序列常常采用
（2）FASTQ格式：表示生物序列及其質(zhì)量的文本格式
（3）ncbi SRA (Sequence Read Archive) ：存儲(chǔ)高通量測序原始數(shù)據(jù)和比對(duì)信息，把FASTQ格式文件壓縮為SRA格式
絕大多數(shù)分析工具不支持SRA，需要使用配套工具包SRA Toolkit先行處理
1. SRA toolkit軟件下載
在官網(wǎng)選擇適合自己的版本下載。
#我選的ubuntu版本，其他一樣，把下載鏈接修改一下就好了
wget http://ftp-trace.ncbi.nlm./sra/sdk/2.10.5/sratoolkit.2.10.5-ubuntu64.tar.gz
用conda install sra-tools失敗，只好用wget方法或者手動(dòng)下載到linux盤符下。把安裝包下載后用tar xzvf 解壓，再配置完PATH就安裝好了。
檢查配置：
prefetch -V
2.用SRAtoolkit下載并處理NCBI數(shù)據(jù)
將 .sra文件轉(zhuǎn)換為 .fstaq.gz文件的工具。用NCBI的SRR數(shù)據(jù)測試一下。
（1）下載
理論上下載東西都可以用wget，但是太慢了。單個(gè)數(shù)據(jù)下載還好，批量下載
prefetch SRRxxxxxxx -O . #-O . 指定到當(dāng)前路徑，否則默認(rèn)路徑難找
一個(gè)數(shù)據(jù)下了好久，大概1個(gè)多小時(shí)。不知道怎么優(yōu)化。
（2）解壓
fastq-dump SRRxxxxxxx.sra #解壓后從sra文件變?yōu)閒astq文件
雙端測序數(shù)據(jù)要加–split-files，否則解壓后兩端的數(shù)據(jù)不會(huì)分開，難以被其他軟件讀取 如果所用分析軟件支持讀取gzip，建議加上–gzip，將解壓后的數(shù)據(jù)用gzip壓縮，避免占用過多空間
fastq-dump --split-files --gzip xxx.sra
（三）測序數(shù)據(jù)質(zhì)控與過濾： fastp
輸出HTML和JSON報(bào)告，前者方便閱讀，后者方便軟件讀取
單端：fastp -i raw.fq -o clean.fq
雙端：fastp -i raw_1.fq -I raw_2.fq -o clean_1.fq -O clean_2.fq
有必要附加的參數(shù)：-l 36 -j xxx.json -h xxx.html
默認(rèn)報(bào)告文件名 fastp.json 和 fastp.html，處理多個(gè)樣本時(shí)極易互相覆蓋，建議改為樣本名稱
fastp參數(shù)設(shè)置
# I/O options 輸入輸出序列文件
-i <單端-輸入文件名>
-o <單端-輸出文件名>
-I <雙端-輸入文件名>
-O <雙端-輸出文件名>
#過濾后的最短序列長度
-l 36 #默認(rèn)15，建議設(shè)為36或40
# reporting options 報(bào)告參數(shù)
-j <the json format report file name >
-h <the html format report file name >
-R "report_title"
二、序列比對(duì)：HISAT2
- 注釋格式介紹
（1）GFF/GTF格式：一般用于基因組和基因注釋
（2）SAM格式：儲(chǔ)存序列比對(duì)到基因組上的信息的文本格式，
（3）BAMS：SAM的基礎(chǔ)上，用二進(jìn)制 (Binary) 編碼，以便壓縮體積。
壓縮率高于gzip，絕大多數(shù)下游分析工具使用
（4）CRAM：在BAM的基礎(chǔ)上，借助參考序列，進(jìn)一步減少空間占用
用SAMtools將SAM轉(zhuǎn)化為BAM或CRAM格式
samtools sort -o xxx.bam xxx.sam
samtools sort -o xxx.cram --reference ref.fa -O cram xxx.sam #加-O指定輸出格式
建立索引以便快速讀取
samtools index xxx.bam
samtools index xxx.cram
- 為什么要比對(duì) (align / map)
locate：測序所得的短序列在基因組的哪個(gè)位置
variant：如果個(gè)別堿基與基因組不一致，是測序錯(cuò)誤還是變異 - 比對(duì)軟件工作過程
根據(jù)基因組序列FASTA和注釋GTF，通過一定的算法編制索引
FASTQ比對(duì)到索引，生成SAM文件
如HISAT 和 Bowtie 基于BWT算法。
1. 用HISAT2建立索引
有注釋：基因組GTF文件Splice Sites和Exons信息，與基因組序列一起用于建立索引
hisat2_extract_splice_sites.py genes.gtf > splicesites.txt
hisat2_extract_exons.py genes.gtf > exons.txt
hisat2-build --ss splicesites.txt --exon exons.txt genome.fa genome
沒注釋：直接用基因組序列建立索引
hisat2-build genome.fa genome
結(jié)果產(chǎn)生索引文件genome（指向.ht結(jié)尾幾個(gè)文件）
2. 比對(duì)
需要用-x指定基因組索引（genome）、-U或者-1、-2輸入FASTQ文件、-S輸出SAM文件，最好還有日志。
hisat2 -x [index location] -U xxx.fq -S xxx.sam --summary-file xxx.align.log --new-summary #單端
hisat2 -x [index location] -1 xxx_1.fq -2 xxx_2.fq -S xxx.sam --summary-file xxx.align.log --new-summary #雙端
比對(duì)結(jié)果解讀
Aligned concordantly：兩端都能合理地比對(duì)上
Aligned discordantly：兩端都比對(duì)上但不合理（位置或方向等不匹配）
unpaired reads：只有一端比對(duì)上
3. 比對(duì)結(jié)果評(píng)估
reads匹配百分比
reads隨機(jī)性分布（reads比對(duì)到基因上的分布均勻說明打斷的隨機(jī)性好）
匹配reads的GC含量和PCR偏好相關(guān)
傳統(tǒng)基于比對(duì)-組裝的方法bam
四、表達(dá)定量和差異表達(dá)分析
（一）表達(dá)水平估計(jì)
在獲得轉(zhuǎn)錄組測序結(jié)果中的轉(zhuǎn)錄本及其功能注釋信息后，就要根據(jù)測序reads比對(duì)到每個(gè)轉(zhuǎn)錄本中的數(shù)目計(jì)算該基因的表達(dá)水平，從而進(jìn)行后續(xù)的分析。
- 表達(dá)定量方法的兩大陣營
（1）Alignment-based
傳統(tǒng)方法，以BAM文件輸入
比對(duì)到基因組：Cufflinks, StringTie，結(jié)果易受測序片段長度影響
比對(duì)到轉(zhuǎn)錄組：eXpress, Salmon，多做一次比對(duì)耗時(shí)偏多
（2）Alignment-free
以FASTQ文件輸入
quasi-mapping ≠ alignment，速度快
結(jié)果較不易受測序片段長度影響
代表工具：kallisto, Sailfish, Salmon
拓展文獻(xiàn)：Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie, and Ballgown
Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis
1.Salmon (Quasi) 流程
salmon也可用于bam輸入，此處以fasta輸入為例：
（1）用salmon index（支持讀取gzip）建立索引
salmon index -t transcripts.fa -i transcripts_index
#可以是fa或fa.gz文件，建立的索引文件為transcripts_index
（2）定量salmon quant分雙端和單端，輸入索引文件transcripts_index，輸出結(jié)果文件夾transcripts_quant
#雙端測序
salmon quant -i transcripts_index -l <LIBTYPE> -1 reads1.fq -2 reads2.fq --validateMappings -o transcripts_quant
#單端測序
salmon quant -i transcripts_index -l <LIBTYPE> -r reads.fq --validateMappings -o transcripts_quant
### --validateMappings 是官方推薦必加參數(shù)，先用敏感策略發(fā)現(xiàn)潛在mapping位點(diǎn)，然后打分并驗(yàn)證，提高準(zhǔn)確度
注意LIBTYPE參數(shù)（1-3位字母）設(shè)置（讓mapping rate正常）：
（3）結(jié)果解讀
輸出文件夾中的quant.sf，是一個(gè)TSV文件。
#EffectiveLength：計(jì)算得到的有效長度，考慮因素包括片段長度分布和序列特異性偏差等，有些下游分析會(huì)用到
#NumReads ：比對(duì)上的reads數(shù)量估計(jì)值，比對(duì)到多處的reads會(huì)根據(jù)相對(duì)豐度產(chǎn)生小數(shù)
#TPM (Transcripts Per Million)：轉(zhuǎn)錄本的相對(duì)豐度估計(jì)值，可用于下游分析
note：轉(zhuǎn)錄本標(biāo)準(zhǔn)化/相對(duì)定量單位
原始的read counts，處理為FPKM，RPKM，TPM等……
三者區(qū)分？什么時(shí)候使用哪個(gè)指標(biāo)？要看清軟件輸入用的指標(biāo)。
（二）差異表達(dá)分析(鑒定差異基因）
1.差異表達(dá)分析的方法和原理
需要將定量后的結(jié)果（表達(dá)矩陣）作為輸入，設(shè)置好分組信息，再進(jìn)行差異表達(dá)分析。
（1）方法：
基于組裝：Cuffdiff, Ballgown，準(zhǔn)確性不足
基于計(jì)數(shù)：DESeq2, edgeR（limma），前者更準(zhǔn)確，后者支持無重復(fù)樣本
→差異表達(dá)分析拓展
其他：GEO2R（針對(duì)GEO數(shù)據(jù)）
（2）標(biāo)準(zhǔn)化
RNA-Seq分析需要對(duì)基因或轉(zhuǎn)錄本的read counts進(jìn)行normalization，因?yàn)槁湓谝粋€(gè)region內(nèi)的read counts取決于基因長度和測序深度。
→拓展文獻(xiàn)Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-Seq data
2.DESeq2流程
（1）準(zhǔn)備輸入文件
①樣本信息矩陣ColData：sample，condition
設(shè)計(jì)比較矩陣（contrast matrix)告訴差異分析函數(shù)應(yīng)該如何對(duì)哪個(gè)因素進(jìn)行比較，[默認(rèn)首字母靠前的condition為對(duì)照!]
②表達(dá)矩陣countData：gene，sample，counts
如果用Salmon、Sailfish、kallisto 得到表達(dá)矩陣，那么就可以用DESeqDataSetFromTximport() 輸入countData。其他導(dǎo)入方法還有DESeqDataSetFromMatrix()、DESeqDataSet()等
- 準(zhǔn)備表達(dá)矩陣實(shí)例：從salmon到tximport
#導(dǎo)入salmon定量的結(jié)果
files <- file.path(samples$run, "quant.sf") #files是一個(gè)個(gè)quants.sf的路徑，選樣本名run一列
#輸入基因ID-TXNAME對(duì)應(yīng)文件
tx2gene <- read.table(file = "tx2gene.txt", sep = "\t")
#定量化生成表達(dá)矩陣
library(tximport)
txi <- tximport(files, type="salmon", tx2gene=tx2gene)
其中，tx2gene是轉(zhuǎn)錄本與基因的轉(zhuǎn)換關(guān)系，可通過AnnotationHub包獲?。?/p>
ah <- AnnotationHub() #下載數(shù)據(jù)庫
sc <- query(ah, 'Saccharomyces cerevisiae') #查詢物種
sc_tx <- sc[['AH64985']] #選擇ID下載詳細(xì)內(nèi)容
k <- keys(sc_tx, keytype = "GENEID") #以基因ID為鍵名
df <- select(sc_tx, keys=k, keytype = "GENEID",columns = "TXNAME") #調(diào)換順序以符合tximport要求：tx2gene <- df[,2:1]
（2）生成DESeqDataSet對(duì)象（tximport 導(dǎo)入為例）
library(DESeq2)
dds <- DESeqDataSetFromTximport(countData, colData = colData, design = ~ condition)
#condition是數(shù)據(jù)框的因子。design說明要分析的變量
#~在R里面用于構(gòu)建公式對(duì)象，~左邊為因變量，右邊為自變量
（3）DESeq2差異表達(dá)分析
①標(biāo)準(zhǔn)化：DESeq()
包括estimation of size factors（estimateSizeFactors)， estimation of dispersion（estimateDispersons)， Negative Binomial GLM fitting and Wald statistics（nbinomWaldTest）三步
dds <- DESeq(dds) #對(duì)dds矩陣對(duì)rawcount進(jìn)行Normalize，不需事先標(biāo)準(zhǔn)化
res <- results(dds) #生成結(jié)果，一個(gè)DESeqResults對(duì)象
summary(res) #用summary看上調(diào)下調(diào)比例（默認(rèn)KD vs control）、離群值等
# resOrdered <- res[order(res$padj),] #p值排序
②可視化：plotMA()
MA圖：M表示log fold change，衡量基因表達(dá)量變化，上調(diào)還是下調(diào)。A表示每個(gè)基因的count的均值。
plotMA(res, ylim=c(-2,2)) #沒有經(jīng)過 statistical moderation平緩log2 fold changes的情況
library(apeglm)
resLFC <- lfcShrink(dds, coef="condition_WT_vs_KD", type="apeglm") #經(jīng)過lfcShrink 收縮log2 fold change
plotMA(resLFC, ylim=c(-2,2))
③確定閾值，篩選差異表達(dá)基因
一般p-value<0.05是顯著, 顯著性不代表結(jié)果正確，只用于給后續(xù)的富集分析和GSEA提供排序標(biāo)準(zhǔn)和篩選。
假陽性隨檢驗(yàn)次數(shù)增加而增加，通常取p<0.05，1000次檢驗(yàn)可以有50次假陽性 Bonferroni
校正：p值除以檢驗(yàn)次數(shù)，0.05/1000=5×10-5，過于嚴(yán)苛導(dǎo)致大量假陰性 False Discovery Rate，常用
Benjamini-Hochberg 即 BH 校正方法 將一系列的p值按照從大到小排序，然后利用公式計(jì)算每個(gè)p值所對(duì)應(yīng)的FDR值：FDR
= p×(n/i)， p是p值，n是p值個(gè)數(shù)，最大的p值的i值為n，第二大則是n-1，依次至最小為1 將計(jì)算出來的FDR值作為新p值，如果某一個(gè)p值所對(duì)應(yīng)的FDR值大于前一位p值（更大的p值）所對(duì)應(yīng)的FDR值，則放棄公式計(jì)算出來的FDR值，選用與它前一位相同的值，因此會(huì)產(chǎn)生連續(xù)相同F(xiàn)DR值的現(xiàn)象；反之則保留計(jì)算的FDR值
返回p值對(duì)應(yīng)的FDR值
res05 <- results(dds, alpha=0.05) #默認(rèn)FDR小于0.1，現(xiàn)取閾值padj小于0.05。padj就是用BH對(duì)多重試驗(yàn)進(jìn)行矯正
res05
summary(res05)
篩選差異顯著的數(shù)據(jù)后，建立基因-FC列表，用作后續(xù)富集分析：
#提取差異表達(dá)基因集：選取上調(diào)FC>2（即log2FC>1）或下調(diào)<-2的基因
diff_gene_info <- subset(res05, (log2FoldChange > 1 | log2FoldChange < -1))
diff_genes <- row.names(diff_gene_info) #
#提取log2FoldChange信息的列表
diff_gene_table <- as.data.frame(diff_gene_info)
geneList <- diff_gene_table[,2]
#log2FoldChange列表用names備注對(duì)應(yīng)基因名稱，排序
names(geneList) = as.character(row.names(diff_gene_table))
geneList <- sort(geneList, decreasing = TRUE)
如果只提取上調(diào)/下調(diào)，步驟也相同，總之DESeq2用于提取我們所需的基因集。
3.edgeR&limma流程
見文章
五、富集分析
富集分析在之前芯片數(shù)據(jù)分析基因的差異表達(dá)的文章中也有寫到，再貼一遍富集分析介紹。
（一）GO富集分析
1.什么是GO（Gene Ontology）
基因已知的功能信息可以分為細(xì)胞組成 Cellular Component (CC)、分子功能 Molecular Function (MF)、生物過程 Biological Process (BP)三個(gè)域。
每一個(gè)域根據(jù)具體功能不同又分為不同 GO term， 有三種關(guān)系：is a，part of，regulates，通過有向無環(huán)圖連接成網(wǎng)
通過分析一組差異基因在功能的分類關(guān)系，可以找到差異基因在那些GO分類條目富集，尋找不同樣品的差異基因可能和哪些基因功能的改變有關(guān)。
官網(wǎng)有詳細(xì)介紹和GO富集分析在線工具。
2.實(shí)現(xiàn)工具
3.GO富集分析：clusterProfiler包
（1）enrichGO()生成enrichResult對(duì)象
輸入：
①待富集的基因集（如差異分析一步得到的上調(diào)基因）
不難理解這種只用了基因集的富集分析算法屬于過表達(dá)分析（over representation analysis）
②物種基因數(shù)據(jù)庫（OrgDb查詢）
library("clusterProfiler")
library("org.xxx.db") #物種基因數(shù)據(jù)庫
enrichGO_up_BP <- enrichGO(gene = up_genes, OrgDb = "org.Sc.sgd.db", keyType = "ENSEMBL", ont = "BP")
#keyType和比對(duì)GTF一致，ont三選一
（2）富集分析結(jié)果可視化
用enrichplot包實(shí)現(xiàn)條形圖barplot()、散點(diǎn)圖dotplot()、有向無環(huán)圖plotGOgraph()的繪制：
library("enrichplot")
barplot(enrichGO_up_BP, showCategory = 20) #條形圖
dotplot(enrichGO_up_BP, showCategory = 20) #散點(diǎn)圖
plotGOgraph(enrichGO_up_BP) #有向無環(huán)圖，顏色表示顯著性，紅色為最顯著的10個(gè)
ggupset包繪制upset圖對(duì)基因集合可視化：
library("ggupset")
upsetplot(enrichGO_up_BP) #upset plot是高階的venn圖，揭示基因和基因集之間的關(guān)系
對(duì)于表達(dá)水平，可以用heatplot()繪制熱圖：
heatplot(enrichGO_up_BP, foldChange = gene_FC_list) #foldChange是排序后的FC-基因列表
（二）KEGG富集分析
1.什么是KEGG PATHWAY
- Kyoto Encyclopedia of Genes and Genomes （KEGG）京都基因與基因組百科全書
- KEGG PATHWAY： is a collection of manually drawn pathway maps representing our knowledge on the molecular interaction, reaction and relation networks for: ①M(fèi)etabolism， ②Genetic Information Processing ，③Environmental Information Processing ，④Cellular Processes ，⑤Organismal Systems，⑥Human Diseases，⑦Drug Development
2.工具
（1）在線工具
KOBAS、
（2）本地工具
clusterProfiler包
3.KEGG富集分析：clusterProfiler包
還是用這個(gè)包，與GO富集分析類似做法，只不過函數(shù)是enrichKEGG(),organism走（物種縮寫查詢）。
enrichKEGG_up <- enrichKEGG(gene = up_genes, organism = "sce", keyType = 'kegg')
barplot(enrichKEGG_up)
dotplot(enrichKEGG_up)
note：著名的clusterProfiler包可以完成許多類富集分析，有空仔細(xì)研究。 →clusterProfiler包富集分析
參考文獻(xiàn)：
網(wǎng)絡(luò)資料、上機(jī)課課件
|
