edger和deseq2_转录组分析(二)Hisat2+DESeq2/EdgeR
一、序列比对
在2016年的一篇综述A survey of best practices for RNA-seq data analysis,提到目前有三种RNA数据分析的策略。那个时候的工具也主要用的是TopHat,STAR和Bowtie.其中TopHat目前已经被它的作者推荐改用HISAT进行替代。
1. Hisat2教程
1.1 下载安装
#conda直接安装
conda install hisat2
#源码下载安装
wget wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/downloads/hisat2-2.1.0-source.zip
unzip hisat2-2.1.0-source.zip
make
1.2 构建index
直接下载现有的insex或通过Hisat2的方法进行创建
# 其实hisat2-buld在运行的时候也会自己寻找exons和splice_sites,但是先做的目的是为了提高运行效率
extract_exons.py gencode.v26lift37.annotation.sorted.gtf > hg19.exons.gtf &
extract_splice_sites.py gencode.v26lift37.annotation.gtf > hg19.splice_sites.gtf &
# 建立index, 必须选项是基因组所在文件路径和输出的前缀
hisat2-build --ss hg19.splice_sites.gtf --exon hg19.exons.gtf genome/hg19/hg19.fa hg19
1.3正式比对
hisat2基本用法就是hisat2 [options]* -x {-1 -2 | -U } [-S ],基本就是提供index的位置,PE数据或者是SE数据存放位置。然而其他可选参数却是进阶的一大名堂。新手就用默认参数呗。
hisat2 --dta -p 6 --max-intronlen 5000000 -x Oryza_sativa.IRGSP-1.0.genome -1 C1-1_good_1.fq -2 C1-1_good_2.fq -S C1-1.HISAT_aln.sam >hisat2_running.log 2>&1
1.4 Hisat2输出结果
比对之后会输出如下结果,解读一下就是全部数据都是100%的,2.88%的配对数据一次都没有比对,94.20%的数据比是唯一比对,2.92%是多个比对。然后如果不按照顺序来,有4.96%的比对。之后把剩下的部分用单端数据进行比对的话,65.57%数据没比对上,33.23%的数据比对一次,1.20%比对超过一次。零零总总的加起来是98.20%的比对。
20182824 reads; of these:
20182824 (100.00%) were paired; of these:
581893 (2.88%) aligned concordantly 0 times
19011569 (94.20%) aligned concordantly exactly 1 time
589362 (2.92%) aligned concordantly >1 times
----
581893 pairs aligned concordantly 0 times; of these:
28886 (4.96%) aligned discordantly 1 time
----
553007 pairs aligned 0 times concordantly or discordantly; of these:
1106014 mates make up the pairs; of these:
725197 (65.57%) aligned 0 times
367552 (33.23%) aligned exactly 1 time
13265 (1.20%) aligned >1 times
98.20% overall alignment rate
2. SAMtools三板斧
SAM(sequence Alignment/mapping)数据格式是目前高通量测序中存放比对数据的标准格式,当然他可以用于存放未比对的数据。而目前处理SAM格式的工具主要是SAMTools
view: BAM-SAM/SAM-BAM 转换和提取部分比对
sort: 比对排序
merge: 聚合多个排序比对
index: 索引排序比对
faidx: 建立FASTA索引,提取部分序列
tview: 文本格式查看序列
#最常用的三板斧就是格式转换,排序,索引。
samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam
samtools sort SRR35899${i}.bam -o SRR35899${i}_sorted.bam
samtools index SRR35899${i}_sorted.bam
3. BAM/SAM文件格式
SAM文件主要由两个部分构成:
header:标记了该SAM文件的一些基本信息,比如版本、按照什么方式排序的、Reference信息等等。
本体:每行为一个reads,不同列记录了不同的信息,列与列之间通过tab分隔。
每列的含义:
MAPQ值:
表示为mapping的质量值,mapping Quality, It equals -10log10Pr{mapping position is wrong}, rounded to the nearest integer, A value 255 indicates that the mapping quality is not available. 该值的计算方法是mapping的错误率的-10log10值,之后四舍五入得到的整数,如果值为255表示mapping值是不可用的,如果是unmapped read则MAPQ为0,一般在使用bwa mem或bwa aln(bwa 0.7.12-r1039版本)生成的sam文件,第五列为60表示mapping率最高,一般结果是这一列的数值是从0到60,且0和60这两个数字出现次数最多
想把小于2的都丢弃:
samtools view -bSq 2 file.sam > filtered.bam
flag的含义:
1 : 代表这个序列采用的是PE双端测序
2: 代表这个序列和参考序列完全匹配,没有插入缺失
4: 代表这个序列没有mapping到参考序列上
8: 代表这个序列的另一端序列没有比对到参考序列上,比如这条序列是R1,它对应的R2端序列没有比对到参考序列上
16:代表这个序列比对到参考序列的负链上
32 :代表这个序列对应的另一端序列比对到参考序列的负链上
64 : 代表这个序列是R1端序列, read1;
128 : 代表这个序列是R2端序列,read2;
256: 代表这个序列不是主要的比对,一条序列可能比对到参考序列的多个位置,只有一个是首要的比对位置,其他都是次要的
512: 代表这个序列在QC时失败了,被过滤不掉了(# 这个标签不常用)
1024: 代表这个序列是PCR重复序列(#这个标签不常用)
2048: 代表这个序列是补充的比对(#这个标签具体什么意思,没搞清楚,但是不常用)
cigar的含义:
cigar中会包含数字,代表了特定match持续了多少nt;以及不同的字符,代表了不同的match情况。
30S512M216N12S (30nt soft clip -> 512nt exact match -> 216nt skipped region -> 12nt soft clip)
30S (30nt soft clip)
40M (40nt exact match)
其中不同的字符及其含义如下:
参考:
https://www.jianshu.com/p/a584d31418f3
https://www.jianshu.com/p/9c87bba244d8
二、htseq-count的使用
HTSeq作为一款可以处理高通量数据的python包,由Simon Anders, Paul Theodor Pyl, Wolfgang Huber等人携手推出HTSeq — A Python framework to work with high-throughput sequencing data。自发布以来就备受广大分析人员青睐,其提供了许多功能给那些熟悉python的大佬们去自信修改使用,同时也兼顾着给小白们提供了两个可以拿来可用的可执行文件 htseq-count(计数) 和 htseq-qa(质量分析)。
具体参考:
https://www.cnblogs.com/triple-y/p/9338890.html
https://blog.csdn.net/herokoking/article/details/78257714
三、基因差异表达分析
1. DESeq2(DESeq2不支持无生物学重复的数据)
library("DESeq2")
#directory
directory
sampleFiles
#sampleFiles
sampleCondition
sampleTable
#sampleTable
ddsHTSeq
#ddsHTSeq
colData(ddsHTSeq)$condition
dds
res
res
#head(res)
png(file="plotMA.png",type="cairo")
plotMA(dds)
dev.off()
mcols(res,use.names=TRUE)
write.csv(as.data.frame(res),file='deseq2.csv')
png(file="plotDispEsts.png",type="cairo")
plotDispEsts(dds)
dev.off()
#pheatmap
sum(res$padj < 0.05, na.rm=TRUE)
library("pheatmap")
select
nt
log2.norm.counts
df
pdf('heatmap.pdf',width = 6, height = 7)
pheatmap(log2.norm.counts, cluster_rows=TRUE, show_rownames=FALSE,cluster_cols=TRUE, annotation_col=df)
dev.off()
2. EdgeR
count.matrix
ID B.count wt.count
AT1G01010 384 314
AT1G01020 661 861
AT1G01030 48 47
AT1G01040 5608 6206
AT1G01046 77 82
AT1G01050 2889 2357
AT1G01060 6039 6296
AT1G01070 408 240
AT1G01073 0 0
edgeR.R
library(edgeR)
data = read.table('count.matrix', header=T, row.names=1, com='')
col_ordering = c(1,2)
rnaseqMatrix = data[,col_ordering]
rnaseqMatrix = round(rnaseqMatrix)
rnaseqMatrix = rnaseqMatrix[rowSums(cpm(rnaseqMatrix) > 1) >= 2,]
conditions = factor(c(rep('B', 1), rep('wt', 1))) #对样本进行分组
exp_study = DGEList(counts=rnaseqMatrix, group=conditions) #建立基因表达列表,利用DEGList函数, 参数counts为read数矩阵,group为上一步的分组变量
exp_study = calcNormFactors(exp_study) #计算各样本内的标准化因子
exp_study_bcv = exp_study
bcv = 0.01 #估计生物学变异系数
et = exactTest(exp_study_bcv, dispersion = bcv ^ 2)
#gene = decideTestsDGE(et, p.value = 0.05, lfc = 0)
#summary(gene)
#head(res)
tTags = topTags(et,n=NULL)
result_table = tTags$table
result_table = data.frame(sampleA='B', sampleB='wt', result_table)
result_table$logFC = -1 * result_table$logFC
write.csv(result_table,file = 'B-wt.csv')
source('Rna-seq/trinityrnaseq-Trinity-v2.4.0/Analysis/DifferentialExpression/R/rnaseq_plot_funcs.R')
pdf('B-wt.pdf')
plot_MA_and_Volcano(result_table$logCPM, result_table$logFC, result_table$FDR)
dev.off()
edger和deseq2_转录组分析(二)Hisat2+DESeq2/EdgeR相关推荐
- edger多组差异性分析_R语言利用edgeR package进行基因差异表达分析 举例
R语言利用edgeR package进行基因差异表达分析 举例 实验数据: 同一组织,分为两组,control vs treat,每组7例sample.数据第一列为基因名,后14列为对应的count. ...
- 转录组分析_转录组分析 | 使用Stringtie对数据进行下游处理
TCGA | GEO | 文献阅读 | 数据库 | 理论知识 R语言 | Bioconductor | 服务器与Linux 接前文: 转录组分析 | fastqc进行质控与结果解读 转录组分析 | 使 ...
- 跟着Cell学单细胞转录组分析(八):单细胞转录组差异基因分析及多组结果可视化
接着单细胞下游分析: 从Cell学单细胞转录组分析(一):开端!!! 跟着Cell学单细胞转录组分析(二):单细胞转录组测序文件的读入及Seurat对象构建 跟着Cell学单细胞转录组分析(三):单细 ...
- 生信步骤|转录组测序上游分析:hisat2+samtools+stringtie
转录组分析在当下研究功能基因组领域十分常用.相关软件组合种类也十分丰富,本文采用了hisat2+samtools+stringtie策略从转录组数据中挖掘差异表达基因.在这里小编整理了一下此套组合的执 ...
- 转录组分析---Hisat2+StringTie+Ballgown使用
转录组分析---Hisat2+StringTie+Ballgown使用 (2016-10-10 08:14:45) 转载▼ 标签: 生物信息学 转录组 1.Hisat2建立基因组索引: First ...
- 39个工具,120种组合深度评估 (转录组分析工具哪家强)
前言 NGS系列文章包括NGS基础.转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这).ChIP-seq分析 (ChIP-seq基本分析流程).单细胞测序分析 (重磅综述:三万字 ...
- edger多组差异性分析_edgeR差异基因分析的一般过程
基于转录组测序获得的定量表达值,识别差异表达变化的基因或其它非编码RNA分子,实际上方法还是非常多的.但就目前来看,DESeq2和edgeR是出现频率最高的两种方法了. DESeq2已经在上一篇文章中 ...
- 在Windows下完成转录组分析
纯粹就是记录下如何在Window下少写代码完成RNA-seq,并没有详细讲解转录组分析的原理和每一步背后的含义.想深入了解可以翻生信技能树,生信宝典,组学大讲堂等等公众号.本文里的命令都是用power ...
- 本周开课 | 第 17 期高级转录组分析和R数据可视化火热报名中!!!
福利公告:为了响应学员的学习需求,经过易生信培训团队的讨论筹备,现安排<高级转录组分析和R数据可视化>于2022年4月15-17 线上/线下课程 (线上课是通过腾讯会议实时直播线下课,实时 ...
最新文章
- 把eclipse从英文调整为中文
- component是什么接口_阿里高级技术专家:整洁的应用架构“长”什么样?
- boost库学习入门篇
- [机器学习-原理篇]支持向量机(SVM)深入理解
- [HNOI2015]开店(树剖+主席树+标记永久化)
- 2022最新最全升级版【精品工具】用Appuploader发布上传iOS APP上架流程简单快速
- 磁盘分区 如何分出整数
- 解决Debian 11系统缺少无线网卡固件rtl8192cfw.bin
- 激光雷达--C16镭神16线三维激光雷达介绍
- IT4IT的前世今生
- 数据库技术之MySQL高级
- 浅谈大数据平台架构设计
- centos7 下mono安装
- Java并发编程-Exchange
- 伯德图 matlab,matlab画三维伯德图,bode图
- 整车试验数采设备技术方案
- React中使用useState数据异步问题解决方法
- CentOS 7安装MySQL集群-GALERA CLUSTER 4 FOR MYSQL 8 RELEASE
- Windows掉激活或重装后激活失败
- data,bss和rodata段的区别与联系
热门文章
- 数据库前台界面设计html,简单Web界面开发(数据库)
- android9.0来电无法获取处理
- 移动开发 网络流量精简攻略
- 关于stimusoft report 报表一行记录有的高度高,有的高度低,错落的问题解决
- java枚举比较大小写_Java 枚举(enum)的学习
- probuff html5,上分神器!Reno5 Pro自带Buff,游戏手感没谁了
- git合并分支(一看就懂)
- lsmod | grep kvm无输出的解决办法
- 【LeetCode】206. 反转链表
- tars框架 php,tars框架安装