- STAR+RSEM+Deseq2視頻課程鏈接已經(jīng)公開(kāi), https://www.bilibili.com/video/BV1KJ411p7WN?p=8
B站上有大量生信干貨, 學(xué)生信看B站足以 , 別問(wèn)我B站是什么哈 ;
生信技能樹(shù)論壇(http://www.biotrainee.com)
2個(gè)樣品2個(gè)重復(fù)(2*2)
1_1
1_2
2_1
2_2
實(shí)用小工具
SAM flags查詢（sam格式文件）
準(zhǔn)備工作
RNA-seq流程原理
1.準(zhǔn)備工作
Using of software
FastQC Multiqc Trimmomatic STAR RSEM Subread HTSeq kallisto Deseq2 ...（conda可以安裝）
查看軟件是否可用踏拜，及時(shí)更新;準(zhǔn)備好文件粟誓；
~/.bashrc
screen -S han #創(chuàng)建新窗口
# 會(huì)話為dead狀態(tài)。使用screen -wipe命令清除該會(huì)話
# 鑒定是雙端還是單端：計(jì)算fq.gz文件的reads;兩個(gè)gz文件的reads數(shù)目相等，為雙端測(cè)序裹赴；
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值檢測(cè)
2.質(zhì)量控制
質(zhì)量控制
ls *gz |xargs -I [] echo 'nohup fastqc [] &'>fastqc.sh
bash fastqc.sh
Fastqc 堿基質(zhì)量分布圖
橫坐標(biāo)代表每個(gè)每個(gè)堿基的位置窑眯，反映了讀長(zhǎng)信息，比如測(cè)序的讀長(zhǎng)為150bp,橫坐標(biāo)就是1到150湘换；
縱坐標(biāo)代表堿基質(zhì)量值宾舅，
圖中的箱線圖代表在每個(gè)位置上所有堿基的質(zhì)量值分布统阿，
中間的紅線代表的是中位數(shù)
用黃色填充的區(qū)域的上下兩端分別代表上四分位數(shù)和下四分位數(shù)；
箱線圖最上方的短線代表90%筹我，最下方的短線代表10%
藍(lán)色的線代表平均值
背景色從上到在下依次為green, orange, red; 分別代表very good, reasonable, poor;將堿基質(zhì)量分成3個(gè)不同的標(biāo)準(zhǔn)
當(dāng)有一個(gè)位置的10%四分位數(shù)小于10或者中位數(shù)小于25時(shí)會(huì)給出警告扶平；
當(dāng)有一個(gè)位置的10%四分位數(shù)小于5或者中位數(shù)小于20時(shí)會(huì)提示失敗蔬蕊；
如下圖所示：
從上面圖中我們可以看出讀長(zhǎng)為75bp;前15bp左右測(cè)序質(zhì)量非常好结澄；
隨著測(cè)序的進(jìn)行，由于試劑的消耗等原因岸夯，測(cè)序質(zhì)量開(kāi)始逐漸降低麻献，最后的60-75bp質(zhì)量就非常的差；
質(zhì)量過(guò)濾
#以一個(gè)樣品的兩個(gè)重復(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.序列比對(duì)
STAR+RSEM 先比對(duì)再進(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長(zhǎng)度-1（150-1）
序列比對(duì)排序
# STAR align 簡(jiǎn)單版本 樣本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 #將上面基因組比對(duì)的信息轉(zhuǎn)化為轉(zhuǎn)錄本比對(duì)的信息,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 查看文件比對(duì)信息
less **.final.out
4.定量
測(cè)試rsem是否安裝
2.png
#samtools 查看bam 文件;（每個(gè)基因上有多少個(gè)reads 的統(tǒng)計(jì)信息的bam文件;）
#為什么tab文件中后面有三個(gè)文件;
從原始的read信息到比對(duì)文件的過(guò)程;
samtools view sample2_Aligned.sortedByCoord.out.bam |head
轉(zhuǎn)錄組數(shù)據(jù)計(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 通過(guò)什么來(lái)進(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 不比對(duì)源譬，直接定量
#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
矩陣合并
#刪除未檢測(cè)到表達(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)前目錄；
打開(kāi)一個(gè)文件夾中的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)
# 控制條件：因子（對(duì)照2喉祭，實(shí)驗(yàn)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)
#實(shí)際包含3步泛烙，請(qǐng)自行學(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)過(guò)矯正后的表達(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()
#畫(huà)火山圖
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
