RNA-Seq數(shù)據(jù)分析：cutadapt+hisat2+samtools+stringtie+deseq2

參考教程：

數(shù)據(jù)：線蟲轉(zhuǎn)錄組測序PF data（雙端測序）;

流程：

  1. cutadapt去接頭;
  2. hisat2建立索引（indexing）;
  3. hisat2+samtools序列比對（alignment）;

通過index進行比對可加快比對速度。

  1. FastQC對當(dāng)前測序數(shù)據(jù)的質(zhì)量進行評估;
  2. Stringtie對數(shù)據(jù)進行拼接、定量（expression）;

一個基因可能由幾段不相鄰的序列組成。

  1. prepDE.py腳本提取基因表達矩陣;
  2. DESeq2差異分析;
  3. Metascape路徑分析差異基因富集（Pathway and Process Enrichment Analysis）。

cutadapt去接頭（adapter trimming）

已知接頭序列read1_adapt、read2_adapt，需要將read1_adapt 的反向互補序列read1_adapt_REV（在線轉(zhuǎn)換）

cutadapt -a read2_adapt -A read1_adapt_REV -m 20 --pair-filter=both -o out_fq1 -p out_fq2 fq1 fq2
#-m：表示去除接頭后如果read長度小于這個m值就不要了

hisat2建立索引（indexing）

  • hisat2官網(wǎng)
    官網(wǎng)提供了一些物種的index文件，如人類、小鼠
  • 基因注釋：Caenorhabditis_elegans.gtf
  • 參考基因組：Caenorhabditis_elegans.WBcel235.dna.toplevel.fa
#加入-ss, -exon，需要消耗上百G的內(nèi)存，不加也可，具體區(qū)別不太清楚
hisat2_extract_splice_sites.py Caenorhabditis_elegans.gtf > Caenorhabditis_elegans.ss
hisat2_extract_exons.py Caenorhabditis_elegans.gtf > Caenorhabditis_elegans.exon
hisat2-build -p 20 Caenorhabditis_elegans.WBcel235.dna.toplevel.fa --ss Caenorhabditis_elegans.ss --exon Caenorhabditis_elegans.exon Caenorhabditis_elegans

#-p 線程數(shù)
#最后一項為輸出的index文件名

hisat2+samtools序列比對（alignment）

hisat2 -x Caenorhabditis_elegans(index_name) -p 20 -1 N1_R1.fastq.gz -2 N1_R2.fastq.gz -S N1.sam

輸出：

26484487 reads; of these:
  26484487 (100.00%) were paired; of these:
    2355777 (8.89%) aligned concordantly 0 times
    23433490 (88.48%) aligned concordantly exactly 1 time
    695220 (2.63%) aligned concordantly >1 times
    ----
    2355777 pairs aligned concordantly 0 times; of these:
      605991 (25.72%) aligned discordantly 1 time
    ----
    1749786 pairs aligned 0 times concordantly or discordantly; of these:
      3499572 mates make up the pairs; of these:
        2745181 (78.44%) aligned 0 times
        713028 (20.37%) aligned exactly 1 time
        41363 (1.18%) aligned >1 times
94.82% overall alignment rate

輸出sam文件，將其轉(zhuǎn)為bam文件（sam的二進制格式），并對bam文件進行排序。

samtools sort -@ 8 -o *.bam *.sam

FastQC

對測序數(shù)據(jù)（去接頭）的質(zhì)量進行評估。

fastqc N1.bam N2.bam N3.bam T1.bam T2.bam T3.bam

Stringtie對數(shù)據(jù)進行拼接、定量（expression）

stringtie -p 8 -G /……/Caenorhabditis_elegans.gtf -e -B -o N1/transcripts.gtf -A N1/gene_abundances.tsv /……/N1.bam

設(shè)置參數(shù)-B，輸出的結(jié)果為ballgown所需要的格式，需要python腳本prepDE.py，才能得到基因表達矩陣。

使用prepDE.py腳本提取基因表達矩陣

進入ballgown文件夾（由stringtie生成，包含所有樣本），將prepDE.py腳本拷貝至當(dāng)前文件夾，

cp /……/prepDE.py ./

使用python命令直接運行腳本，注意需要使用python2

python prepDE.py

運行結(jié)果中"gene_count_matrix.csv"是基因表達矩陣（DESeq2的輸入文件）。

刪除gene_count_matrix.csv中全是0的行（R語言）：

gene_count <- read.csv("gene_count_matrix.csv",stringsAsFactors = F,header = T)
rownames(gene_count) <- gene_count[,1]
gene_count <-gene_count[,-1]
gene_count<-gene_count[which(rowSums(gene_count) > 0),]

DESeq2差異分析

基于負二項式分布，包含標(biāo)準化（基于負二項式分布）。

參考教程：

# 設(shè)置工作目錄
setwd("C:/Users/DELL/Desktop/gene_count/")
# 讀入基因表達值，設(shè)定行名為gene_id
gene_count <- read.csv("gene_count_matrix.csv",stringsAsFactors = F)
# 對gene_id一列進行拆分，去除重復(fù)名稱
library(stringr)
# 設(shè)置空的"gene_count_1"向量，行數(shù)與上面的測序結(jié)果一致
gene_count_1<-rep(NA,nrow(gene_count))
# 利用for循環(huán)，對gene_count數(shù)據(jù)框中的重復(fù)列進行拆分提取
#將gene_id中|兩側(cè)拆分、提取
for (i in 1:nrow(gene_count)){
  gene_count_1[i] <- unlist(str_split(gene_count[i,1], pattern = "\\|"))[1]
}
# 對原數(shù)據(jù)框中的特定序列重新賦值
gene_count$gene_id <- gene_count_1
# 將第一列作為文件的行名
rownames(gene_count) <- gene_count[,1]
gene_count <-gene_count[,-1]
#保存文件至對應(yīng)目錄
write.csv(gene_count, file = "……/gene_count/gene_count.csv", row.names = TRUE)

#加載DESeq2包
library(DESeq2)
#DESeq2需要三種矩陣，分別為countData表達矩陣,colData樣品信息矩陣及design差異表達矩陣
#countData為表達矩陣即gene_count
#colData為樣品信息矩陣即coldata
#design為差異表達矩陣即批次和條件（對照、處理）等
#設(shè)置condition樣品組別、重復(fù)數(shù)
condition<- factor(c(rep("N", 3), rep("T", 3)), levels = c("N","T"))
condition
#設(shè)置樣品信息矩陣colData
colData <- data.frame(row.names = colnames(gene_count), condition)
colData

#在R里面用于構(gòu)建公式對象，~左邊為因變量，右邊為自變量。
#標(biāo)準流程：dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design= ~ batch + condition) 
#countData為表達矩陣即countdata
#colData為樣品信息矩陣即coldata
#design為差異表達矩陣即批次和條件（對照、處理）等
#對dds進行標(biāo)準流程構(gòu)建
dds <- DESeqDataSetFromMatrix(gene_count, colData, design = ~condition)
#過濾低豐度數(shù)據(jù)
dds <- dds[rowSums(counts(dds)) > 1, ]
#或者在構(gòu)建dds之前加上gene_count <- gene_count[apply(gene_count, 1, sum) > 1 , ] 
#對原始dds進行normalize
dds <- DESeq(dds)
#顯示dds信息
dds

# 對差異分析結(jié)果進行保存 -------------------------------------------------------------

#使用DESeq2包中的results()函數(shù)，提取差異分析的結(jié)果
#Usage:results(object, contrast, name, .....）
#將提取的差異分析結(jié)果定義為變量"res" 
#contrast: 定義誰和誰比較，處理組在前，對照組在后
#提取分析結(jié)果并保存為res
res = results(dds, contrast=c("condition","T","N"))
#對結(jié)果res利用order()函數(shù)按pvalue值進行排序
#order()函數(shù)先對數(shù)值排序，然后返回排序后各數(shù)值的索引，常用用法：V[order(V)]或者df[order(df$variable),]
#對res進行排序
res = res[order(res$pvalue),]
#顯示res結(jié)果前6行信息
head(res)
#對res矩陣進行總結(jié)，利用summary命令統(tǒng)計顯示一共多少個genes上調(diào)和下調(diào)
summary(res)
#將差異分析的所有結(jié)果進行輸出保存
write.csv(res, file="……/all_different_genes.csv")

#利用table函數(shù)統(tǒng)計顯著差異基因的數(shù)目
table(res$padj<0.05&abs(res$log2FoldChange>1))
#對具有顯著性差異的結(jié)果進行過濾、提取
#獲取padj小于0.05，表達倍數(shù)取以2為對數(shù)后的絕對值大于1
#使用subset()函數(shù)過濾需要的結(jié)果至新的變量significant_different_genes_group中
#Usage:subset(x, ...)，其中x為objects，...為篩選參數(shù)或條件
#對group中數(shù)據(jù)進行過濾、提取
significant_different_genes <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
#使用dim函數(shù)查看該結(jié)果的維度、規(guī)模
dim(significant_different_genes)
#顯示結(jié)果的前6行信息
head(significant_different_genes)
#對顯著差異基因進行輸出保存
write.csv(significant_different_genes, file = "……/significant_different_genes.csv")

Metascape路徑分析差異基因富集（Pathway and Process Enrichment Analysis）

上調(diào)基因（up）：log2FoldChange>0

step2：選擇物種；
step3：express analysis（快速分析），custom analysis（自定義分析）。

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容