轉(zhuǎn)錄組分析——STAR+RSEM+Deseq2

B站上有大量生信干貨, 學(xué)生信看B站足以 , 別問我B站是什么哈 ;
生信技能樹論壇(http://www.biotrainee.com)

2個樣品2個重復(fù)(2*2)
1_1
1_2
2_1
2_2

實用小工具

SAM flags查詢（sam格式文件）


準(zhǔn)備工作



RNA-seq流程原理
1.準(zhǔn)備工作
Using of software

FastQC Multiqc Trimmomatic STAR RSEM Subread HTSeq kallisto Deseq2 ...（conda可以安裝）

查看軟件是否可用，及時更新;準(zhǔn)備好文件；
 ~/.bashrc
screen -S han   #創(chuàng)建新窗口
# 會話為dead狀態(tài)。使用screen -wipe命令清除該會話

# 鑒定是雙端還是單端：計算fq.gz文件的reads;兩個gz文件的reads數(shù)目相等，為雙端測序；
zcat -f gKO-01.fq.gz |wc -l

# 參考基因組和注釋下載；(示例而已)
[ftp://ftp.ensembl.org/pub/release-95/](ftp://ftp.ensembl.org/pub/release-95/)

# Genome file
[ftp://ftp.ensembl.org/pub/release-95/fasta/mus_musculus/dna/](ftp://ftp.ensembl.org/pub/release-95/fasta/mus_musculus/dna/)[Mus_musculus.GRCm38.dna.toplevel.fa.gz](ftp://ftp.ensembl.org/pub/release-95/fasta/mus_musculus/dna/Mus_musculus.GRCm38.dna.toplevel.fa.gz)

# GTF file
[ftp://ftp.ensembl.org/pub/release-95/gtf/mus_musculus/](ftp://ftp.ensembl.org/pub/release-95/gtf/mus_musculus/)[Mus_musculus.GRCm38.95.gtf.gz](ftp://ftp.ensembl.org/pub/release-95/gtf/mus_musculus/Mus_musculus.GRCm38.95.gtf.gz)

# 解壓參考基因組和注釋文件信息
gunzip *.gz

# 解壓縮命令
tar -xzvf file.tar.gz //解壓tar.gz

# count md5 value,，查看文件是否完整 
cd 01raw_data
ls
md5sum *gz>md5.txt
cat md5.txt
ls
md5sum -c md5.txt
MD5值檢測
2.質(zhì)量控制
質(zhì)量控制
ls *gz |xargs -I [] echo 'nohup fastqc [] &'>fastqc.sh
bash fastqc.sh
Fastqc 堿基質(zhì)量分布圖

橫坐標(biāo)代表每個每個堿基的位置，反映了讀長信息，比如測序的讀長為150bp,橫坐標(biāo)就是1到150；
縱坐標(biāo)代表堿基質(zhì)量值，
圖中的箱線圖代表在每個位置上所有堿基的質(zhì)量值分布，
中間的紅線代表的是中位數(shù)
用黃色填充的區(qū)域的上下兩端分別代表上四分位數(shù)和下四分位數(shù)；
箱線圖最上方的短線代表90%，最下方的短線代表10%
藍(lán)色的線代表平均值
背景色從上到在下依次為green, orange, red; 分別代表very good, reasonable, poor;將堿基質(zhì)量分成3個不同的標(biāo)準(zhǔn)
當(dāng)有一個位置的10%四分位數(shù)小于10或者中位數(shù)小于25時會給出警告；
當(dāng)有一個位置的10%四分位數(shù)小于5或者中位數(shù)小于20時會提示失?。?br> 如下圖所示：



從上面圖中我們可以看出讀長為75bp;前15bp左右測序質(zhì)量非常好；
隨著測序的進(jìn)行，由于試劑的消耗等原因，測序質(zhì)量開始逐漸降低，最后的60-75bp質(zhì)量就非常的差；

質(zhì)量過濾

#以一個樣品的兩個重復(fù) 為例；
sample1_1_R1.fastq.gz sample1_1_R2.fastq.gz \
trimmomatic PE -threads 4 \
sample1_1_R1.fastq.gz sample1_1_R2.fastq.gz \
../02clean_data/sample1_1_paired_clean_R1.fastq.gz \
../02clean_data/sample1_1_unpair_clean_R1.fastq.gz \
../02clean_data/sample1_1_paired_clean_R2.fastq.gz \
../02clean_data/sample1_1_unpair_clean_R2.fastq.gz \
ILLUMINACLIP:/home/miniconda3/share/trimmomatic-0.36-5/adapters/TruSeq3-PE-2.fa:2:30:10:1:true \
LEADING:3 TRAILING:3 \
SLIDINGWINDOW:4:20 MINLEN:50 TOPHRED33

sample1_2_R1.fastq.gz sample1_2_R2.fastq.gz \
trimmomatic PE -threads 10 \
sample1_2_R1.fastq.gz sample1_2_R2.fastq.gz \
../02clean_data/sample1_2_paired_clean_R1.fastq.gz \
../02clean_data/sample1_2_unpair_clean_R1.fastq.gz \
../02clean_data/sample1_2_paired_clean_R2.fastq.gz \
../02clean_data/sample1_2_unpair_clean_R2.fastq.gz \
ILLUMINACLIP:/home/miniconda3/share/trimmomatic-0.36-5/adapters/TruSeq3-PE-2.fa:2:30:10:1:true \
LEADING:3 TRAILING:3 \
SLIDINGWINDOW:4:20 MINLEN:50 TOPHRED33
3.序列比對
STAR+RSEM 先比對再進(jìn)行定量


# build index 建立索引;
STAR --runThreadN 10 --runMode genomeGenerate \
--genomeDir arab_STAR_genome \
--genomeFastaFiles 00ref/TAIR10_Chr.all.fasta \
--sjdbGTFfile 00ref/Araport11_GFF3_genes_transposons.201606.gtf \
--genomeSAindexNbases 12 \
--sjdbOverhang 149  # 一般為reads長度-1（150-1）

序列比對排序
# STAR align 簡單版本 樣本R1
STAR --runThreadN 10 --genomeDir arab_STAR_genome \          #索引目錄
--readFilesCommand zcat \       #執(zhí)行解壓縮：
--readFilesIn 02clean_data/sample1_1_paired_clean_R1.fastq.gz \
02clean_data/sample1_1_paired_clean_R2.fastq.gz \          #執(zhí)行輸入文件
--outFileNamePrefix 03align_out/sample1_1 \       #輸出文件加前綴；
--outSAMtype BAM SortedByCoordinate \          #指定輸出bam 文件并排序；
--outBAMsortingThreadN 5 \                              #設(shè)置排序中的線程數(shù)；
--quantMode TranscriptomeSAM GeneCounts #將上面基因組比對的信息轉(zhuǎn)化為轉(zhuǎn)錄本比對的信息,count
# --quantMode TranscriptomeSAM will output alignments translated into  tran-script coordinates,為了使用RSEM 進(jìn)行定量分析做準(zhǔn)備；

 # STAR align 復(fù)雜版本 樣本R2
STAR --runThreadN 12 --genomeDir arab_STAR_genome \
--readFilesCommand zcat \
--readFilesIn 02clean_data/sample2_1_paired_clean_R1.fastq.gz \
02clean_data/sample2_1_paired_clean_R2.fastq.gz \
--outFileNamePrefix 03align_out/sample2_1 \
--outSAMtype BAM SortedByCoordinate \
--outBAMsortingThreadN 5 \
--quantMode TranscriptomeSAM GeneCounts

#less 查看文件比對信息
less **.final.out

4.定量

測試rsem是否安裝

2.png

#samtools 查看bam 文件;（每個基因上有多少個reads 的統(tǒng)計信息的bam文件;）
#為什么tab文件中后面有三個文件;
從原始的read信息到比對文件的過程;
samtools view sample2_Aligned.sortedByCoord.out.bam |head

轉(zhuǎn)錄組數(shù)據(jù)計算表達(dá)定量;
# rsem prepare reference  從參考基因組中提取原始準(zhǔn)備文件
rsem-prepare-reference --gtf 00ref/Araport11_GFF3_genes_transposons.201606.gtf \
00ref/TAIR10_Chr.all.fasta \
arab_RSEM/arab_rsem
RSEM建立索引
#sample
rsem-calculate-expression --paired-end --no-bam-output \
--alignments -p 15 \
-q 03align_out/sample1_1Aligned.toTranscriptome.out.bam \
arab_RSEM/arab_rsem \
04rsem_out/sample1_1rsem
RSEM進(jìn)行定量
基因組，轉(zhuǎn)錄本表達(dá)定量結(jié)果
查看RSEM定量結(jié)果
基因表達(dá)定量方案選擇
featureCount定量
featureCount
###featurecount
## programname -p(雙端) -a (基因組注釋文件) -o (輸出文件) -T(使用的線程數(shù)) -t -g 通過什么來進(jìn)行定量，輸出什么結(jié)果   所有要定量的文件;
featureCounts -p -a ../00ref/Araport11_GFF3_genes_transposons.201606.gtf -o out_counts.txt -T 6 -t exon -g gene_id sample*Aligned.sortedByCoord.out.bam
featureCounts -p -a ../00ref/genes.gtf -o Z-13-1_out_counts.txt -T 10 -t exon -g *L1*dm6_Aligned.sortedByCoord.out.bam
featureCounts -p -a ../00ref/genes.gtf -o out_counts.txt -T 20 -t exon -g gene_id

featureCounts -p -a ../00ref/genes.gtf -o out_counts.txt -T 6 -t exon -g gene_id zgq-*_accepted_hits.bam

方案

kallisto 不比對，直接定量

#kallisto  先構(gòu)建索引 index;   /arab_kallisto目錄下
kallisto index -i arab_kallisto ../arab_RSEM/arab_rsem.transcripts.fa

#kallisto進(jìn)行定量
kallisto quant -i Dro_kallisto/Dro_kallisto -o 05kallisto_out/sample2 02clean_data/R1.gz 02clean_data/R2.gz
kallisto quant -i Dro_kallisto/Dro_kallisto -o 05kallisto_out/Z-1 ../../Data_cy/20181018-embryo-RIP-seq\ data/Z-13-1_L1_310310.R1.clean.fastq.gz ../../Data_cy/20181018-embryo-RIP-seq\ data/Z-13-1_L1_310310.R2.clean.fastq.gz 
#abundance.tsv    列：gene   counts   TPM
kallisto建立索引

kallisto定量結(jié)果

5.差異基因分析：


R install

install.packages("ggrepel")
install.packages("ggplot2")


###  差異表達(dá)
mkdir 06deseq_out

#構(gòu)建表達(dá)矩陣
rsem-generate-data-matrix *rsem.genes.results > output.matrix
矩陣合并
#刪除未檢測到表達(dá)的基因
awk 'BEGIN{printf "geneid\ta1\ta2\tb1\tb2\n"}{if($2+$3+$4+$5>0)print $0}' output.matrix > deseq2_input.txt
###查看文件的行數(shù)
wc -l output.matrix 
wc -l deseq2_input.txt 

mv deseq2_input.txt ../06deseq_out/
mv deseq2.r ../06deseq_out/

abundance_estimates_to_matrix.pl
run_DE_analysis.pl
差異基因結(jié)果
差異基因作圖

MAplot

火山圖

R腳本

### 設(shè)置當(dāng)前目錄；
打開一個文件夾中的R腳本，session
或 輸入路徑；setwd("~/00")



# read.table ,文件有header,第一行為row.name;
input_data <- read.table("deseq2_input.txt", header=TRUE, row.names=1)

# 取整,函數(shù)round
input_data <-round(input_data,digits = 0)

# 準(zhǔn)備文件
# as.matrix 將輸入文件轉(zhuǎn)換為表達(dá)矩陣；
input_data <- as.matrix(input_data)

# 控制條件：因子（對照2，實驗2）
condition <- factor(c(rep("ct1", 2), rep("exp", 2)))

# input_data根據(jù)控制條件構(gòu)建data.frame
coldata <- data.frame(row.names=colnames(input_data), condition)

library(DESeq2)
# build deseq input matrix 構(gòu)建輸入矩陣
#countData作為矩陣的input_data；colData Data.frame格式；控制條件design;
dds <- DESeqDataSetFromMatrix(countData=input_data,colData=coldata, design=~condition)

# DESeq2 進(jìn)行差異分析
dds <- DESeq(dds)

#實際包含3步，請自行學(xué)習(xí);
# 提取結(jié)果
#dds dataset格式轉(zhuǎn)換為DESeq2 中result數(shù)據(jù)格式，矯正值默認(rèn)0.1
res <- results(dds,alpha=0.05)
#查看res(DESeqresults格式),可以看到上下調(diào)基因
summary(res)
上下調(diào)基因
#res(resultset)按照校正后的P值排序
res <- res[order(res$padj), ]
res

# 將進(jìn)過矯正后的表達(dá)量結(jié)果加進(jìn)去；
resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)),by="row.names",sort=FALSE)
names(resdata)[1] <- "Gene"

#查看(resdata)
# output result 輸出結(jié)果

write.table(resdata,file="diffexpr-results.txt",sep = "\t",quote = F, row.names = F)

#可視化展示
plotMA(res)
maplot <- function (res, thresh=0.05, labelsig=TRUE,...){
with(res,plot(baseMean, log2FoldChange, pch=20, cex=.5, log="x", ...))
with(subset(res, padj<thresh), points(baseMean, log2FoldChange, col="red", pch=20,cex=1.5))

}
png("diffexpr-malot.png",1500,1000,pointsize=20)
maplot(resdata, main="MA Plot")
dev.off()

#畫火山圖
install.packages("ggrepel")
library(ggplot2)
library(ggrepel)

resdata$significant <- as.factor(resdata$padj<0.05 & abs(resdata$log2FoldChange) > 1)
ggplot(data=resdata, aes(x=log2FoldChange, y =-log10(padj),color =significant)) +
geom_point() +
ylim(0, 8)+
scale_color_manual(values =c("black","red"))+
geom_hline(yintercept = log10(0.05),lty=4,lwd=0.6,alpha=0.8)+
geom_vline(xintercept = c(1,-1),lty=4,lwd=0.6,alpha=0.8)+
theme_bw()+
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"))+
labs(title="Volcanoplot_biotrainee (by sunyi)", x="log2 (fold change)",y="-log10 (padj)")+

theme(plot.title = element_text(hjust = 0.5))+
geom_text_repel(data=subset(resdata, -log10(padj) > 6), aes(label=Gene),col="black",alpha =0.8)

差異基因的提取
#查看符合條件的差異基因
awk '{if($3>1 && $7<0.05)print $0}' diffexpr-results.txt |head
#查看差異基因有多少行
awk '{if($3>1 && $7<0.05)print $0}' diffexpr-results.txt |wc -l
#提取某幾列
awk '{if($3>1 && $7<0.05)print $0}' diffexpr-results.txt |cut -f 1,2,4,7 |head
#上調(diào)基因的提取;
awk '{if($3>1 && $7<0.05)print $0}' diffexpr-results.txt |cut -f 1,2,4,7 > up.gene.txt
#下調(diào)基因的提取;
awk '{if($3<-1 && $7<0.05)print $0}' diffexpr-results.txt |cut -f 1,2,4,7 > down.gene.txt

根據(jù)第1列是Geneid，第7,8，9，10列是counts數(shù)，用awk提取出geneID和counts。
awk -F '\t' '{print $1,$7,$8,$9,$10}' OFS='\t' out_counts.txt >subread_matrix.out
最后編輯于 ：
?著作權(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)容