转录组分析之RNA-seq实战

参考：https://github.com/shenmengyuan/RNA_seq_Biotrainee

基础知识部分

Fastqc 碱基质量分布图

横坐标代表每个每个碱基的位置，反映了读长信息，比如测序的读长为150bp,横坐标就是1到150；

纵坐标代表碱基质量值，

图中的箱线图代表在每个位置上所有碱基的质量值分布，

中间的红线代表的是中位数

用黄色填充的区域的上下两端分别代表上四分位数和下四分位数；

箱线图最上方的短线代表90%，最下方的短线代表10%

蓝色的线代表平均值

背景色从上到在下依次为green, orange, red; 分别代表very good, reasonable, poor;将碱基质量分成3个不同的标准

当有一个位置的10%四分位数小于10或者中位数小于25时会给出警告；

当有一个位置的10%四分位数小于5或者中位数小于20时会提示失败；

如下图所示：

 从上面图中我们可以看出读长为75bp;前15bp左右测序质量非常好；

随着测序的进行，由于试剂的消耗等原因，测序质量开始逐渐降低，最后的60-75bp质量就非常的差；

# count md5 value,，查看文件是否完整

cd 01raw_data

ls

md5sum *gz>md5.txt

cat md5.txt

ls

md5sum -c md5.txt

质量控制：

ls *gz |xargs -I [] echo 'nohup fastqc [] &'>fastqc.sh

bash fastqc.sh

质量过滤：

trimmomatic PE -threads 4 \

sample1_R1.fastq.gz sample1_R2.fastq.gz \

../02clean_data/sample1_paired_clean_R1.fastq.gz \

../02clean_data/sample1_unpair_clean_R1.fastq.gz \

../02clean_data/sample1_paired_clean_R2.fastq.gz \

../02clean_data/sample1_unpair_clean_R2.fastq.gz \

ILLUMINACLIP:/home/sy/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

STAR+RSEM 先比对再进行定量

# build index 建立索引;

STAR --runThreadN 6 --runMode genomeGenerate \

--genomeDir arab_STAR_genome \

--genomeFastaFiles 00ref/TAIR10_Chr.all.fasta \

--sjdbGTFfile 00ref/Araport11_GFF3_genes_transposons.201606.gtf \

--sjdbOverhang 149 # 一般为reads长度-1（150-1）

序列比对排序

# STAR align 简单版本 样本R1

STAR --runThreadN 5 --genomeDir arab_STAR_genome \

--readFilesCommand zcat \

--readFilesIn 02clean_data/sample1_paired_clean_R1.fastq.gz \

02clean_data/sample1_paired_clean_R2.fastq.gz \

--outFileNamePrefix 03align_out/sample1 \

--outSAMtype BAM SortedByCoordinate \

--outBAMsortingThreadN 5 \

--quantMode TranscriptomeSAM GeneCounts

# STAR align 复杂版本 样本R2

STAR --runThreadN 5 --genomeDir arab_STAR_genome \

--readFilesCommand zcat \

--readFilesIn 02clean_data/sample2_paired_clean_R1.fastq.gz \

02clean_data/sample2_paired_clean_R2.fastq.gz \

--outFileNamePrefix 03align_out/sample2 \

--outSAMtype BAM SortedByCoordinate \

--outBAMsortingThreadN 5 \

--quantMode TranscriptomeSAM GeneCounts

# rsem prepare reference  从参考基因组中提取原始准备文件

rsem-prepare-reference --gtf 00ref/Araport11_GFF3_genes_transposons.201606.gtf \

00ref/TAIR10_Chr.all.fasta \

arab_RSEM/arab_rsem

rsem-calculate-expression --paired-end --no-bam-output \

--alignments -p 5 \

-q 03align_out/sample2_Aligned.toTranscriptome.out.bam \

arab_RSEM/arab_rsem \

04rsem_out/sample2_rsem

基因组，转录本表达定量结果

#htseq-count 未进行实验

htseq-count -r pos -m union -f bam -s no \

-q 03align_out/sample2Aligned.sortedByCoord.out.bam

05htseq_out/sample2.htseq.out

kallisto 不比对，直接定量

#kallisto  先构建索引 index;

kallisto index -i arab_kallisto ../arab_RSEM/arab_rsem.transcripts.fa

差异表达分析

###featurecount

## programname -p(双端) -a (基因组注释文件) -o (输出文件) -T(使用的线程数) -t -g 通过什么来进行定量，输出什么结果  所有要定量的文件;

feartureCounts -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

###  差异表达

mkdir 06deseq_out

#构建表达矩阵

rsem-generate-date-matrix *_rsem.genes.results > output.matrix

#删除未检测到表达的基因

awk 'BEGIN{printf "geneid\ta1\ta2\tb1\tb2\n"}{if($2+$3+$4+$5>0)print $0}'

output.matrix > deseq2_input.txt

###

abundance_estimates_to_matrix.pl

run_DE_analysis.pl

Rstudio

MA plot 图形

差异基因的提取

#查看符合条件的差异基因

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

#上调基因的提取;

awk '{if($3>1 && $7<0.05)print $0}' diffexpr-results.txt |cut -f 1,2,4,7 > up.gene.txt

#下调基因的提取;

awk '{if($3<-1 && $7<0.05)print $0}' diffexpr-results.txt |cut -f 1,2,4,7 > down.gene.txt

附录

R脚本

---

#

input_data <- read.table("deseq2_input.txt", header=TRUE, row.names=1)

#取整

input_data <-round(input_data,digits = 0)

# r准备文件

input_data <- as.matrix(input_data)

condition <- factor(c(rep("ct1", 2), rep("exp", 2)))

coldata <- data.frame(row.names=colnames(input_data), condition)

library(DESeq2)

# build deseq input matrix 构建输入矩阵

dds <- DESeqDataSetFromMatrix(countData=input_data,colData=coldata, design=~condition)

#  DESeq2 进行差异分析

dds <- DESeq(dds)

#实际包含3步，请自行学习;

# 提取结果

res <- results(dds,alpha=0.05)

summary(res)

res <- res[order(res$padj), ]

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 输出结果

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

}

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)