【生信技能树】转录组测序数据分析

视频地址转录组测序数据分析

2020-12-2完成视频的学习及笔记记录；后面学习完下一个视频后，会有一波对笔记的补充

学习小结：这是第二次看视频，第一次就大概浏览了一下，这一次非常认真地做了笔记，收货还是非常大的，了解了整个测序数据的分析过程，从原始的测序数据到后面的一些简单的分析有了更深的理解，对之前分析处理数据的一些疑惑有了解答，算是基本掌握了整个RNA-seq的分析流呈，上游chip和count是没有问题了，下游的差异分析，可视化部分出一些简单的热图、火山图、相关线性图、柱状也基本能搞定；但是其实要学的还有很多，很明显的感觉到问题还存在很多，好好学习过linux和R以后还得再把这个视频过一遍，希望下次看的时候能有更多的收获。

P1：生物学背景知识

幕布

1、一些常识训练 2、转录组背景知识获得：①收集RNA-seq技术综述；②阅读相关RNA-seq的文献，公司的结题报告 3、了解RNA-seq的实验环节：实验设计环节、RNA的提取及质量控制、cDNA的合成、文库构建 4、了解RNA-seq的应用：①蛋白质基因编码结构②新型蛋白质编码基因③基因表达的量化和比较④表达数量性状基因座⑤单细胞RNA-seq⑥融合基因⑦基因变异⑧长的非编码RNA⑨非编码小RNA⑩扩增产物测序 5、根据测序策略及各个分组是否有重复来实战；SE50：无重复；PE150：有重复 6、了解一些实战导读 7、数据处理的流程：质控、比对、计数、差异分析 8、可视化：IGV、ggplot2+ggpubr包

P2：转录组背景知识

1、Linux背景知识 ①系统认知：马哥linux运维 ②去可视化： ③文本处理：主要命令--awk/grep/sed/paste/cat/diff/wc/vi；处理fastq、fasta、bam、sam、gff、gtf、bed、MAF等格式 ④软件安装： 分为三类--第一是二进制可执行软件，直接下载软件包解压后既可以全路径下调用；第二就是所有的语言代码：perl、R、python、java、matlab、ruby、C；第三是各种软件安装器：apt-get、yum、bioconda、cpan、cran、brew、pip、conda（conda贼好用）

2、基础知识介绍掌握 ①编程基础（终身学习）： linux持续学习，马哥视频+练习题 生信人的linux考试 生信分析人员如何系统入门Linux(2019更新版) R语言持续学习，视频+练习题 生信人的20个R语言练习题 生信分析人员如何系统入门R(2019更新版) python或者perl脚本任选一个 要求：给出视频答案并录制视频讲解自己的解题思路 ②安装软件（越多越好） ③生信基础知识掌握（终身学习）生物信息学基础知识100讲 1）生物芯片及测序数据的分类，原理和历史；主要的测序平台；主要的芯片平台 2）三大国际数据中心的了解：NCBI、ENSEMBL、UCSC（可以在谷歌浏览器输入UCSC database filetype：pdf 就能检索到其下数据库的一些介绍） 3）数据格式的整理和熟记：fastq、fasta、sam、bam、vcf、gff、gtf、bed、MAF~~ 4）参考基因组的熟悉及其基因组注释新文件下载及摸索 5）从基因开始理解生物信息学 6）组学技术应用的第一篇文章极易最新综述文章收集整理 7）各个组学数据分析的结题报告的阅读及整理 8）数据库的收集整理，包括

P3：转录组破冰之旅

1、高通量测序方案的选择 转录调控研究：转录组出、表达谱测序、Small RNA测序、CircRNA测序、LncRNA测序、全长转录组测序、甲基化测序 微生物组学研究：环境微生物多样性检测、宏基因组de novo测序、宏转录组测序 基因组学研究：全基因组de novo测序、简化基因组测序、基因组重测序、外显子组测序、扩增子测序

2、实验设计 sample RNA→RNA fragments→ cDNA fragments→reads

3、一些思考 生物学重复比测序深度重要的多，对结果影响更大 转录组分析生物学重复n≥30（经费充足）；经费紧张n≥6；重复少的话假阴性率越高，筛选到的差异表达基因越少。；（一般很少有测到6个以上的。。。。。。） 如果转录组分析重复不足（n＜6）时，需要我门用更严格的分析方法如：DESeq、edgeR、sleuth；

4、核心流程 比对--计数--差异分析 Tophat-Cufflink-Cuffdiff Subread-featureCounts-DESeq2 STAR-RSEM-EBSeq Bowtie-eXpress-edgeR kallisto-sleuth HiSAT-StringTie-Ballgown

对于应用者来说没有必要把每一个软件搞得那么清楚，好好地使用软件得到想要的结果才是王道 没有绝对最优解，统计知识很重要

P4：转录组文献解读

19年：RNA-seq这十年（3万字长文综述) 根据文献中的方法和步骤来不断优化、丰富自己的数据处理

1、数据处理的流程 主要代码都是在： 原创10000+生信教程大神给你的RNA实战视频演练 ①数据资源下载，参考基因组及参考转录组，基因注释文件； 参考基因到ucsc下载？注释文件在ensembl下载？ 注释文件里罗列了基因的转录本、对应的外显子 参考基因组及注释文件下载--浏览器搜索：hg19 ftp ucsc/ncbi/ensembl --选择hg38.fa.gz 右键点击复制下载链接，挂在服务器上下载：wget、url

less 命令可以查看hg19.fa的文件的具体内容，每个基因参考文件·

②质控，需要fastqc及multiqc，trim-galore；或者很多其它软件 trimmomatic, cutadapt, qualimap等等 ③比对： 很多软件：star, hisat2,bowtie2,tophat,bwa,subread，优先选择STAR或者hisat2 ④计数 软件也很多：featureCounts,htseq, bedtools ,deeptools salmon ，优先选择featureCounts ⑤normalization 归一化（统计学+数学知识背景） Total Count (TC)；Upper Quartile (UQ)；Median (Med)；the DESeq normalization implemented in the DESeq Bioconductor package；Trimmed Mean of M values (TMM) implemented in the edgeR Bioconductor package；Quantile (Q)；the Reads Per Kilobase per Million mapped reads (RPKM) normalization, FPKM,TPM ⑥差异分析等 主要是导入R里面使用R包DEseq2或者edgeR分析，也可以选择limma的voom算法 可以参考我GitHub的示例

P5-6：软件安装-上

Conda：目前最强的非root软件安装管理器 如何安装conda： miniconda是conda的一个轻量版本

####下载miniconda
wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-py38_4.9.2-Linux-x86_64.sh 
####安装miniconda
bash Anaconda3-2020.11-Linux-x86_64.sh 
####将miniconda3保存到环境命令中
在 .bashrc 里最后面添加一行命令
export PATH="/home/feng/miniconda3/bin:$PATH"

####安装一个rna的环境
miniconda3/bin/conda create -n rna python=2
conda create -n rna python=2 #创建名为rna的软件安装环境
conda info --envs #查看当前conda环境
source activate rna #激活conda的rna环境

####安装软件前可以搜索一下conda里是否存在这个软件
conda search softwarename 
###安装相关生信软件
conda install -y sra-tools
conda install -y trimmomatic
conda install -y cutadapt multiqc 
conda install -y trim-galore
conda install -y star hisat2 bowtie2
conda install -y subread tophat htseq bedtools deeptools
conda install -y salmon
source deactivate #注销当前的rna环境

P7：数据下载

1、去GEO上获得id清单，然后用prefetch下载

image.png
另一种挂载后台的方法；除了nohup；还可以直接把命令括起来；但是有多少个下载任务就会分配多少个进程，如果要去kill的话，需要一个一个去kill；可以写一个循环，一个一个去kill

下载完一定要检查仔细，看文件大小是否匹配

ps -ef |grep prefe|awk  '{print $2}'|while read id ;do kill $id ;done    ###把进程号全部选中并做一个循环杀死命令

2、sra2fastq

ls sra* |while read id ;do (nohup fastq-dump --gzip --split-3 -o ./ $id &);done       

P8：qc1

检查fastq文件 用zless，因为我们上一步导出的文件是压缩的文件（--gzip），所以需要zless查看

可以看到fastq是每四行一段序列；但是一行一行看太费劲了。

image.png

fastqc -t 10 SRR*  ##-t是指用几个线程
QC完后用multiqc 将所有的qc报告整合在一起看
multiqc ./

P9 qc-2

trim_galore

conda activate rna
cd ../2_fq/
ls *_1.fastq.gz >1
ls *_2.fastq.gz >2
paste 1 2 > config
dir='../4_clean/'
nohup cat config| while read id
do
{
     arr=($id)
     fq1=${arr[0]}    
     fq2=${arr[1]}
     trim_galore -q 25  --phred33 --stringency 3 --paired -o ../4_clean $fq1 $fq2 &    
     joblist=($(jobs -p))
     while ((${#joblist[*]} >= 8 ))
     do
     {
          sleep 1
          joblist=($(jobs -p))
     }
     done
     sleep 10
}
done >../4_clean/nohup.log 2>&1 & 
conda deactivate rna

##jm老师
zcat SRRname |grep TCTTCCTTG      ##查找srr文件里的这段序列 

phred33是测序公司的一种质量值体系，还有一种是phred64，现在大多都是phred33为主。质量体系phred33和phred64的解答 --length 50 是指测序长度去掉左右两边后的长度＜50就舍弃调，这个值是根据bp长度为150来说的；但是对于上面截图出现的length=63时，需要将50改为35，不然会舍弃太多的序列 --stringency 3 多长的序列算接头，需要去除 adatater ；常为3 --paired 双端测序数据

image.png

P10、alignment-1

1、构建索引

jm视频里的代码如下

image.png

--STAR runMode genomeGenerate \    ##构建模式
--genomeDir  ~/reference/index/star/hg19 \   ##输出位置
--genomeFastaFiles ~/reference/index/star/hg19/hg19.fa \   ##通过fastafiles文件来产生索引
--sjdbGTFfile ~/reference/gtf/geneome/hg19/hg19.fa \  ##需要基因转换文件
--sjdbOverhang 149 \  ##测序是150，所以选149
--runThreadN 4  ##给4个线程

各个软件索引的构建的代码都能在网上找到，后续实战的时候会把代码贴上
2、比对

trim过后的文件

image.png

创建比对后的文件夹

image.png

ls *gz|while read id;do(echo {$id%%.*});done     ##%%是把点.去掉

image.png

ls *gz |while read id;do (zcat $id|head -1000 > $(basename $id ".gz"));done  #### 创建文件名字去除.gz的文件夹

①hisat比对
hisat -p -x -1 -2 -S (线程数；索引；输入文件1；输入文件2；输出文件)
②subjunc比对
Subjunc -T -i -r -R -o(线程数；索引；输入文件1；输入文件2；输出文件)
③bowtie2比对
bowtie2 -p -x -1 -2 -S(线程数；索引；输入文件1；输入文件2；输出文件) 比对率下降，不能跨内含子比对
④bwa
bwa mem -t 5 -M 1 2 >

P11、alignment-2

1、sam文件里面有什么
①每条染色体的长度
②会记录输出sam文件的代码，可以偷师
③每条比对序列的质量；
30M33S：比对上了30个碱基，33个碱基没比对上
下图是bwa的sam文件

image.png
下图是hisat的sam里的一段：46M4640N17M：跨越了4640个内含子
image.png
2、找一个sam文件格式解读的文章 博客--SAM格式详解简书SAM/BAM 格式文件内容解析
3、bam文件怎么读--二进制文件
用less直接看的话是一对乱码；这里需要用到samtools工具

samtool view erniufd.bam |less -S       ####读取后用less的sort命令，更方便查看
###samtool tview ####类似于IGV，可以查看整条染色体的位置；进入tview后，用/输入你想看的位置，

P12、alignment-3

1、比对软件之间的结果差别
hisat大一点，bwa和bowtie2差不多；subjuncs明显小于其它，大概1/3，这是因为subjuncs跑出来的其实是一个bam文件（二进制）数据量会小很多，但是老师命名成了sam文件，所以看起来差别这么大
但是文件名错了一点也不重要，如果更改了会影响到后面对数据的批量处理。

image.png
image.png

ls *.sam|head ##输出前十个sam文件名
ls *.sam|xargs head ##输出每个sam文件里的前十行

2、批量转换bam文件

nohup ls *.sam |while read id
do
samtools view -b -@ 10 -o ../7_bam/$(basename ${id} ".sam").bam ${id}
done > ../7_bam/nohup.log 2>&1 & 

3、查看比对成功率
overall alignment rate
4、对bam文件qc

ls *.bam | while read id
do
    samtools flagstat -@ 10 ${id} > ../9_bamqc/$(basename ${id} ".bam").stat
done

5、5个比对软件的bam文件qc比较（flagstat）
同样可以用multiqc出一个报告

P13、表达矩阵探索

1、count

conda activate rna
cd ~/GEO/8_sortbam/
gtf="/home/stu3/data/hg38.knownGene.gtf"
featureCounts -T 10 -p -t exon -g gene_id -a $gtf -o ~/GEO/10_count/all.id.txt *.bam
##-p:双端测序
##-t exon 根据外显子来count
##-g gene_id count后取基因名

R语言处理

##去除gene名的点号
library(stringer)
a2$id <- str_split(a2$id,'\\.'simplify=T)[,1]
tmp<-merge(a1,a2,by='id')
##画相关线性图
png('tmp.png')
plot(tmp(,2:3))
dev.off

image.png

P14：DEG-1

健明老师的代码
数据处理R代码

##去重复
duplicated
##去缺失值
DEG = na.omit(DEG)
##转化为数字型矩阵
exprSet<-as.data.frame(lapply(exprSet,as.numeric),row.names = rownames(aeg1))   
##数字归一化
exprSet<-as.data.frame(lapply(exprSet,scale),row.names = rownames(aeg1)) 

1、三种edgeR、DESeq2、Limma差异分析包的代码

##DESeq2
##需要注意的是：group_list是分组信息；exprSet是表达矩阵
##建议直接导入表达矩阵加分组信息，这样代码可以一键运行到底
#exprSet--表达矩阵；group_list--分组信息
suppressMessages(library(DESeq2))   ##加载包
#group_list=c('untrt','trt' ,'untrt','trt' ,'untrt','trt','untrt','trt')   ##设置分组信息
(colData <- data.frame(row.names=colnames(exprSet),     
                       group_list=group_list) )
dds <- DESeqDataSetFromMatrix(countData = exprSet,          ##dds需要coungtData、colData、design这三个输入数据
                              colData = colData,
                              design = ~ group_list)
dds <- DESeq(dds)
res <- results(dds, 
               contrast=c("group_list","trt","untrt"))     ##group_list里trt比untrt
resOrdered <- res[order(res$padj),]     ##给p值排序
head(resOrdered)
DEG =as.data.frame(resOrdered)
DEG = na.omit(DEG)

##edgeR     这个差异分析讲的不是很详细，有些地方读不懂，得先学一下R，再回来弄懂它，不过代码是能跑下来的
library(edgeR)
#group_list=c('untrt','trt' ,'untrt','trt' ,'untrt','trt','untrt','trt')   ##设置分组信息；exprSet<-分析的矩阵
d <- DGEList(counts=exprSet,group=factor(group_list))
keep <- rowSums(cpm(d)>1) >= 2    ##这里及后面的代码不是很懂，老师也没讲很清楚，得去看看R语言的视屏再回来看看
table(keep)
d <- d[keep, , keep.lib.sizes=FALSE]           
d$samples$lib.size <- colSums(d$counts)
d <- calcNormFactors(d)
d$samples
design <- model.matrix(~0+factor(group_list))
dge=d
rownames(design)<-colnames(dge)
colnames(design)<-levels(factor(group_list))
dge <- estimateGLMCommonDisp(dge,design)
dge <- estimateGLMTrendedDisp(dge, design)
dge <- estimateGLMTagwiseDisp(dge, design)
fit <- glmFit(dge, design)
# https://www.biostars.org/p/110861/
lrt <- glmLRT(fit, contrast=c(0,1))    ##从design里面看，0代表trt，1代表utrt
nrDEG=topTags(lrt, n=nrow(dge))
nrDEG=as.data.frame(nrDEG)
head(nrDEG)
edgeR_DEG =nrDEG

##Limma    这个包老师讲的也不是很详细，学过R以后再回来看看
suppressMessages(library(limma))
#group_list=c('untrt','trt' ,'untrt','trt' ,'untrt','trt','untrt','trt')   ##设置分组信息；exprSet<-分析的矩阵
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
dge <- DGEList(counts=exprSet)
dge <- calcNormFactors(dge)
logCPM <- cpm(dge, log=TRUE, prior.count=3)
v <- voom(dge,design,plot=TRUE, normalize="quantile")
fit <- lmFit(v, design)
group_list
cont.matrix=makeContrasts(contrasts=c('untrt-trt'),levels = design)   ####untrt比上trt
fit2=contrasts.fit(fit,cont.matrix)
fit2=eBayes(fit2)
tempOutput = topTable(fit2, coef='untrt-trt', n=Inf)
DEG_limma_voom = na.omit(tempOutput)

####
save(DEG_limma_voom,DEseq_DEG,edgeR_DEG,       ####将三个差异分析包的结果保存下来，后续可以做一个对比
     dds,exprSet,group_list,
     file = 'DEG_results.Rdata')

差异分析----建立好三个矩阵：表达矩阵、分组矩阵、比较矩阵
列名和行名并不算在行数或列数里面
需要看results的说明书才能弄懂为什么是"trt"比上"untrt"
2、相关线性图----代码可以成功运行，里面的一些细节后续再来完善，尽量把每一句话读懂，这个可能得去学习一下画图代码

getwd()
setwd("/Users/xiaofashi/Desktop/shengxingstudy/DEG /")
if(F){
  library(airway)
  data(airway)
  exprSet=assay(airway)
  group_list=colData(airway)[,3]
  group_list=relevel(group_list,ref = 'trt')
  save(exprSet,group_list,file = 'airway_exprSet.Rdata')
}
####相关性图
###导入一个带表达矩阵及分组信息的数据
rm(list = ls())
options(stringsAsFactors = F)
load(file = 'airway_exprSet.Rdata')
if(T){
  colnames(exprSet)
  pheatmap::pheatmap(cor(exprSet))
  group_list
  tmp=data.frame(g=group_list)
  rownames(tmp)=colnames(exprSet)
  # 组内的样本的相似性应该是要高于组间的！
  pheatmap::pheatmap(cor(exprSet),annotation_col = tmp)
  dim(exprSet)
  exprSet=exprSet[apply(exprSet,1, function(x) sum(x>1) > 5),]
  dim(exprSet)
  
  exprSet=log(edgeR::cpm(exprSet)+1)
  dim(exprSet)
  exprSet=exprSet[names(sort(apply(exprSet, 1,mad),decreasing = T)[1:500]),]
  dim(exprSet)
  M=cor(log2(exprSet+1))
  tmp=data.frame(g=group_list)
  rownames(tmp)=colnames(M)
  pheatmap::pheatmap(M,annotation_col = tmp)
  pheatmap::pheatmap(M,annotation_col = tmp,filename = 'cor.pdf')
  
  #library(pheatmap)
  #pheatmap(scale(cor(log2(exprSet+1))))
}

image.png

3、热图

nrDEG=DEG
library(pheatmap)
choose_gene=head(rownames(nrDEG),50)
choose_matrix=exprSet[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))     ###神奇的操作，加完以后图瞬间变好看了！！！##其实就是归一化--Z-score
pheatmap(choose_matrix,filename='DEG_top100_heatmap.png')

image.png

4、火山图

colnames(DEG)
need_DEG<-DEG
logFC_cutoff <- with(need_DEG,mean(abs(log2FoldChange)) + 2*sd(abs( log2FoldChange)) )
# logFC_cutoff=1
need_DEG$change = as.factor(ifelse(need_DEG$pvalue < 0.05 & abs(need_DEG$log2FoldChange) > logFC_cutoff,
                                   ifelse(need_DEG$log2FoldChange > logFC_cutoff ,'UP','DOWN'),'NOT'))  ####将数据分为三组
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),    
                    '\nThe number of up gene is ',nrow(need_DEG[need_DEG$change =='UP',]) ,
                    '\nThe number of down gene is ',nrow(need_DEG[need_DEG$change =='DOWN',]))
library(ggplot2)
g = ggplot(data=need_DEG,
           aes(x=log2FoldChange, y=-log10(pvalue),
               color=change)) +
  geom_point(alpha=0.4, size=1.75) +
  theme_set(theme_set(theme_bw(base_size=20)))+
  xlab("log2 fold change") + ylab("-log10 p-value") +
  ggtitle( this_tile ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
  scale_colour_manual(values = c('blue','black','red')) ## corresponding to the levels(res$change)
print(g)
ggsave(g,filename ='volcano.png')

image.png

5、一些柱状图

rld <- rlogTransformation(dds)
exprMatrix_rlog=assay(rld)
x=apply(exprMatrix_rlog,1,mean)
y=apply(exprMatrix_rlog,1,mad)
plot(x,y)
png("DESeq2_RAWvsNORM.png",height = 800,width = 800)
par(cex = 0.7)
n.sample=ncol(exprSet)
if(n.sample>40) par(cex = 0.5)
cols <- rainbow(n.sample*1.2)
par(mfrow=c(2,2))
boxplot(exprSet, col = cols,main="expression value",las=2)
boxplot(exprMatrix_rlog, col = cols,main="expression value",las=2)
hist(as.matrix(exprSet))
hist(exprMatrix_rlog)
dev.off()
}

image.png

P15：DEG-2

1、查看三个差异分析之间的差别 ####视频上因为时间问题，没讲完，留了个作业，后面给它解决掉

rm(list = ls())
options(stringsAsFactors = F)
load(file = 'DEG_results.Rdata')
a1=data.frame(gnee=rownames(DEG_limma_voom),limma=DEG_limma_voom$logFC)
a2=data.frame(gnee=rownames(DEseq_DEG),limma=DEseq_DEG$log2FoldChange)
a3=data.frame(gnee=rownames(edgeRDEG),limma=edgeR_DEG$logFC)

发现a3的logFC值过大，明显不对

2^27次方是什么概念，就是表达相差上千万倍

image.png