最新更新见 介绍本文参考 bioconductor 中RNA-Seq workflow: gene-level exploratory analysis and differential expression并对其根据需要进行了增减。
更多细节还请参考 试验数据数据来源
Himes BE, Jiang X, Wagner P, Hu R, Wang Q, Klanderman B, Whitaker RM, Duan Q, Lasky-Su J, Nikolos C, Jester W, Johnson M, Panettieri R Jr, Tantisira KG, Weiss ST, Lu Q. “RNA-Seq Transcriptome Profiling Identifies CRISPLD2 as a Glucocorticoid Responsive Gene that Modulates Cytokine Function in Airway Smooth Muscle Cells.” PLoS One. 2014 Jun 13;9(6):e99625. PMID: . GEO: GSE52778.在这个RNA-Seq试验中，采用了4种呼吸道平滑肌肉细胞（airway smooth muscle cells），每种细胞均有 地塞米松治疗、非治疗两类。共计8个样本，储存在 airway 包中。原始数据的处理高通量测序数据常采用 FASTQ 格式来保 存所测的碱基读段和质量分数。如图 所示,FASTQ 格式以测序读段为单位存 储,每条读段占 4 行,其中第一行和的第三行由文件识别标志和读段名(ID)组成(第一行以“@”开头而第三行以“+”开头;第三行中 ID 可以省略,但“+”不能省 略),第二行为碱基序列,第四行为各碱基所对应的测序质量分数序列。采用 tophat/bowtie2 将原始数据fastq映射到基因组序列，得到bam文件； 此处我们采用airway 自带的bam文件。加载 airway 包, 并利用自带bam文件[code]library("airway")
dir &- system.file("extdata", package="airway", mustWork=TRUE)
list.files(dir) # 文件列表[code]##
[1] "GSE52778_series_matrix.txt"
"Homo_sapiens.GRCh37.75_subset.gtf"
[3] "sample_table.csv"
"SraRunInfo_SRP033351.csv"
[5] "SRR1039508_subset.bam"
"SRR1039508_subset.bam.bai"
[7] "SRR1039509_subset.bam"
"SRR1039512_subset.bam"
[9] "SRR1039513_subset.bam"
"SRR1039516_subset.bam"
## [11] "SRR1039517_subset.bam"
"SRR1039520_subset.bam"
## [13] "SRR1039521_subset.bam"[code]csvfile &- file.path(dir,"sample_table.csv")
(sampleTable &- read.csv(csvfile,row.names=1)) # 获取样本信息[code]##
SampleName
Run avgLength Experiment
## SRR1039508 GSM1275862
N61311 untrt untrt SRR1039508
SRX384345 SRS508568 SAMN
## SRR1039509 GSM1275863
trt untrt SRR1039509
SRX384346 SRS508567 SAMN
## SRR1039512 GSM2611 untrt untrt SRR1039512
SRX384349 SRS508571 SAMN
## SRR1039513 GSM2611
trt untrt SRR1039513
SRX384350 SRS508572 SAMN
## SRR1039516 GSM0611 untrt untrt SRR1039516
SRX384353 SRS508575 SAMN
## SRR1039517 GSM0611
trt untrt SRR1039517
SRX384354 SRS508576 SAMN
## SRR1039520 GSM1011 untrt untrt SRR1039520
SRX384357 SRS508579 SAMN
## SRR1039521 GSM1011
trt untrt SRR1039521
SRX384358 SRS508580 SAMN[code]filenames &- file.path(dir, paste0(sampleTable$Run, "_subset.bam")) # 提取bam文件获取bam数据[code]library("Rsamtools")
filenames[code]## [1] "/Library/Frameworks/R.framework/Versions/library/airway/extdata/SRR1039508_subset.bam"
## [2] "/Library/Frameworks/R.framework/Versions/library/airway/extdata/SRR1039509_subset.bam"
## [3] "/Library/Frameworks/R.framework/Versions/library/airway/extdata/SRR1039512_subset.bam"
## [4] "/Library/Frameworks/R.framework/Versions/library/airway/extdata/SRR1039513_subset.bam"
## [5] "/Library/Frameworks/R.framework/Versions/library/airway/extdata/SRR1039516_subset.bam"
## [6] "/Library/Frameworks/R.framework/Versions/library/airway/extdata/SRR1039517_subset.bam"
## [7] "/Library/Frameworks/R.framework/Versions/library/airway/extdata/SRR1039520_subset.bam"
## [8] "/Library/Frameworks/R.framework/Versions/library/airway/extdata/SRR1039521_subset.bam"[code]bamfiles &- BamFileList(filenames, yieldSize=2000000) #
将bam文件放入列表，yieldSize 表示每次被读取的记录数
seqinfo(bamfiles) # 序列的基本信息[code]## Seqinfo object with 84 sequences from an unspecified genome:
seqlengths isCircular genome
&NA&[code]seqlevels(bamfiles) # 所有的染色体名称， GLxxxxxx.x表示genomic contig[code]##
## [15] "22"
## [22] "9"
"GL" "GL" "GL"
## [29] "GL" "GL" "GL" "GL" "GL" "GL" "GL"
## [36] "GL" "GL" "GL" "GL" "GL" "GL" "GL"
## [43] "GL" "GL" "GL" "GL" "GL" "GL" "GL"
## [50] "GL" "GL" "GL" "GL" "GL" "GL" "GL"
## [57] "GL" "GL" "GL" "GL" "GL" "GL" "GL"
## [64] "GL" "GL" "GL" "GL" "GL" "GL" "GL"
## [71] "GL" "GL" "GL" "GL" "GL" "GL" "GL"
## [78] "GL" "GL" "GL" "GL" "GL" "GL" "GL"导入基因组特征（注释）eg. 外显子的染色体位置， 基因的起始、终止位点[code]library("GenomicFeatures")
gtffile &- file.path(dir,"Homo_sapiens.GRCh37.75_subset.gtf")
(txdb &- makeTxDbFromGFF(gtffile, format="gtf"))[code]## TxDb object:
## # Db type: TxDb
## # Supporting package: GenomicFeatures
## # Data source: /Library/Frameworks/R.framework/Versions/library/airway/extdata/Homo_sapiens.GRCh37.75_subset.gtf
## # Organism: NA
## # miRBase build ID: NA
## # Genome: NA
## # transcript_nrow: 65
## # exon_nrow: 279
## # cds_nrow: 158
## # Db created by: GenomicFeatures package from Bioconductor
## # Creation time:
17:40:06 +0800 (Wed, 17 Jun 2015)
## # GenomicFeatures version at creation time: 1.20.1
## # RSQLite version at creation time: 1.0.0
## # DBSCHEMAVERSION: 1.1[code](genes &- exonsBy(txdb, by="gene"))[code]## GRangesList object of length 20:
## GRanges object with 18 ranges and 2 metadata columns:
ranges strand
| &integer&
&character&
1 [087705]
1 [090307]
1 [090939]
1 [094963]
1 [097868]
1 [107176]
1 [107176]
1 [107280]
1 [107284]
1 [107290]
## &19 more elements&
## -------
## seqinfo: 1 sequence from a no seqlengths[code]seqlevels(genes) # 染色体的名字[code]## [1] "1"染色体的名称 seqlevels(bamfiles) 与 seqlevels(genes) 应该保持一致， 特别留意是否含有 “chr”, 要么都有”chr”, 要么都没有。
注意，这里采用了一个基因组注释文件的子集， 完整的信息可以从
获取基因计数[code]library("GenomicAlignments")
se &- summarizeOverlaps(features = genes, reads = bamfiles,
mode = "Union",
# 读段覆盖的模式
singleEnd=FALSE, #双末端 not 单末端
ignore.strand=TRUE,# True 表示忽略±链的限制
fragments=TRUE
) # 只应用于双末端测序,true表示非成对的对端应该被计数
class(se) # 得到 SummarizedExperiment 数据，可用于后续计算[code]## [1] "SummarizedExperiment"
## attr(,"package")
## [1] "GenomicRanges"上图显示的是SummarizedExperiment类(以及他的子类DESeqDataSet)的布局, 粉红色 assay(se) 表示实际的数据， 每行为一个基因，每列为一个样本； colData 表示样本的具体信息，随后我们会对它进行填充；rowRanges 表示每一个基因的信息。具体如下[code]se[code]## class: SummarizedExperiment
## dim: 20 8
## exptData(0):
## assays(1): counts
## rownames(20): ENSG ENSG ... ENSG ENSG
## rowRanges metadata column names(0):
## colnames(8): SRR1039508_subset.bam SRR1039509_subset.bam ... SRR1039520_subset.bam
SRR1039521_subset.bam
## colData names(0):[code]head(assay(se))[code]##
SRR1039508_subset.bam SRR1039509_subset.bam SRR1039512_subset.bam
SRR1039513_subset.bam SRR1039516_subset.bam SRR1039517_subset.bam
SRR1039520_subset.bam SRR1039521_subset.bam
1590[code]colSums(assay(se))[code]## SRR1039508_subset.bam SRR1039509_subset.bam SRR1039512_subset.bam SRR1039513_subset.bam
## SRR1039516_subset.bam SRR1039517_subset.bam SRR1039520_subset.bam SRR1039521_subset.bam
9168[code]colData(se)[code]## DataFrame with 8 rows and 0 columns[code]rowRanges(se)[code]## GRangesList object of length 20:
## GRanges object with 18 ranges and 2 metadata columns:
ranges strand
| &integer&
&character&
1 [087705]
1 [090307]
1 [090939]
1 [094963]
1 [097868]
1 [107176]
1 [107176]
1 [107280]
1 [107284]
1 [107290]
## &19 more elements&
## -------
## seqinfo: 1 sequence from a no seqlengths[code]str(metadata(rowRanges(se)))[code]## List of 1
$ genomeInfo:List of 14
..$ Db type
: chr "TxDb"
..$ Supporting package
: chr "GenomicFeatures"
..$ Data source
: chr "/Library/Frameworks/R.framework/Versions/library/airway/extdata/Homo_sapiens.GRCh37.75_subset.gtf"
..$ Organism
..$ miRBase build ID
..$ Genome
..$ transcript_nrow
: chr "65"
..$ exon_nrow
: chr "279"
..$ cds_nrow
: chr "158"
..$ Db created by
: chr "GenomicFeatures package from Bioconductor"
..$ Creation time
: chr " 17:40:06 +0800 (Wed, 17 Jun 2015)"
..$ GenomicFeatures version at creation time: chr "1.20.1"
..$ RSQLite version at creation time
: chr "1.0.0"
..$ DBSCHEMAVERSION
: chr "1.1"[code](colData(se) &- DataFrame(sampleTable)) # 填充样本的具体信息，方便后续分组，寻找差异基因[code]## DataFrame with 8 rows and 9 columns
SampleName
Run avgLength Experiment
&factor& &factor& &factor& &factor&
&factor& &integer&
## SRR1039508 GSM1275862
untrt SRR1039508
SRX384345 SRS508568
## SRR1039509 GSM1275863
untrt SRR1039509
SRX384346 SRS508567
## SRR1039512 GSM1275866
untrt SRR1039512
SRX384349 SRS508571
## SRR1039513 GSM1275867
untrt SRR1039513
SRX384350 SRS508572
## SRR1039516 GSM1275870
untrt SRR1039516
SRX384353 SRS508575
## SRR1039517 GSM1275871
untrt SRR1039517
SRX384354 SRS508576
## SRR1039520 GSM1275874
untrt SRR1039520
SRX384357 SRS508579
## SRR1039521 GSM1275875
untrt SRR1039521
SRX384358 SRS508580
## SRR1039508 SAMN
## SRR1039509 SAMN
## SRR1039512 SAMN
## SRR1039513 SAMN
## SRR1039516 SAMN
## SRR1039517 SAMN
## SRR1039520 SAMN
## SRR1039521 SAMN
注意：此处得到的数据需要采用EDSeq2包进行差异分析，所以不对数据进行标准化，切记。差异表达基因分析我们采用 DESeq2 包进行，差异表达基因的分析[code]# 此步采用 airway 包自带的se数据进行后续操作，可以忽略。如果没有进行上面的步骤也可以直接采用下面的数据进行后续操作。
data("airway")
se &- airway[code]library("DESeq2")
dds &- DESeqDataSet(se, design = ~ cell + dex) # design 参数为 formula，此处为cell和dex两个因素，~ cell + dex表示我们想控制cell研究dex的影响。采用DESeqDataSetFromMatrix函数从matrix中获取数据[code]countdata &- assay(se)
# 可以根据自己的需要填充自己的数据（matrix格式），这里以assay(se)为例
class(countdata)
head(countdata)
coldata &- colData(se)
(ddsMat &- DESeqDataSetFromMatrix(countData = countdata,
colData = coldata,
design = ~ cell + dex))[code]dds$dex &- relevel(dds$dex, "untrt") # 将 untrt 定义为dex因素的第一水平，随后的foldchange 将采用 trt/untrt
dds &- DESeq(dds)
(res &- results(dds)) # 得到结果，可以根据padj来挑选合适的差异表达基因，log2FoldChange来确定基因上调还是下调，pvalue的校正采用了Benjamini-Hochberg方法，具体见 ?p.adjust[code]## log2 fold change (MAP): dex trt vs untrt
## Wald test p-value: dex trt vs untrt
## DataFrame with 20 rows and 6 columns
baseMean log2FoldChange
-0...2289484
-1...6580603
## ENSG 34.9593217
-0...8663899
0.[code]mcols(res, use.names=TRUE)[code]## DataFrame with 6 rows and 2 columns
description
&character&
&character&
## baseMean
intermediate mean of normalized counts for all samples
## log2FoldChange
log2 fold change (MAP): dex trt vs untrt
standard error: dex trt vs untrt
Wald statistic: dex trt vs untrt
Wald test p-value: dex trt vs untrt
BH adjusted p-values[code]summary(res)[code]##
## out of 16 with nonzero total read count
## adjusted p-value & 0.1
## LFC & 0 (up)
## LFC & 0 (down)
## outliers [1]
## low counts [2]
## (mean count & 12.1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results保存数据[code]write.csv(res, file = '/your/path/')sessionInfo[code]sessionInfo()[code]## R version 3.2.0 ()
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X 10.10.3 (Yosemite)
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## attached base packages:
## [1] parallel
grDevices utils
## other attached packages:
[1] GenomicAlignments_1.4.1
GenomicFeatures_1.20.1
AnnotationDbi_1.30.1
[4] Biobase_2.28.0
knitr_1.10.5
BiocStyle_1.6.0
[7] Rsamtools_1.20.4
Biostrings_2.36.1
XVector_0.8.0
## [10] airway_0.102.0
DESeq2_1.8.1
RcppArmadillo_0.5.200.1.0
## [13] Rcpp_0.11.6
GenomicRanges_1.20.5
GenomeInfoDb_1.4.0
## [16] IRanges_2.2.4
S4Vectors_0.6.0
BiocGenerics_0.14.0
## [19] readr_0.1.1
sqldf_0.4-10
RSQLite_1.0.0
## [22] DBI_0.3.1
gsubfn_0.6-6
proto_0.3-10
## [25] dplyr_0.4.1
plyr_1.8.3
## loaded via a namespace (and not attached):
[1] splines_3.2.0
Formula_1.2-1
assertthat_0.1
latticeExtra_0.6-26
[5] yaml_2.1.13
lattice_0.20-31
chron_2.3-45
digest_0.6.8
[9] RColorBrewer_1.1-2
colorspace_1.2-6
htmltools_0.2.6
XML_3.98-1.2
## [13] biomaRt_2.24.0
genefilter_1.50.0
zlibbioc_1.14.0
xtable_1.7-4
## [17] snow_0.3-13
scales_0.2.5
BiocParallel_1.2.3
annotate_1.46.0
## [21] ggplot2_1.0.1
nnet_7.3-9
survival_2.38-1
magrittr_1.5
## [25] evaluate_0.7
MASS_7.3-40
foreign_0.8-63
tools_3.2.0
## [29] formatR_1.2
stringr_1.0.0.9000
munsell_0.4.2
locfit_1.5-9.1
## [33] cluster_2.0.1
lambda.r_1.1.7
futile.logger_1.4.1
grid_3.2.0
## [37] RCurl_1.95-4.6
bitops_1.0-6
tcltk_3.2.0
rmarkdown_0.7
## [41] gtable_0.1.2
reshape2_1.4.1
gridExtra_0.9.1
rtracklayer_1.28.4
## [45] Hmisc_3.16-0
futile.options_1.0.0 stringi_0.4-1
geneplotter_1.46.0
## [49] rpart_4.1-9
acepack_1.3-3.3[code]library(knitr)
knit('/Users/lipidong/baiduyun/work/RFile/MarkDown/funSet.Rmd', output = '~/learn/blog/_posts/-RNA-Seq.md')fold change
芯片筛选显示,模型组高表达miRNA106个,低表达miRNA 63个, 差异倍数 ( fold change )在1.8～7.0。
基于6个网页-
基因4 × 44 K 芯片进行检测。采用GeneSpring11.0 统计软件，以P 倍数变化（fold change，FC）≥ 2 筛选出差异表达基因，并对部分差异表达基因行RT-PCR 验证。 结果 成功建立SCI 动物模
基于1个网页-
采用倍数变化分析法
MMP-3、MMP-9在创伤性深静脉血栓形成中作用的实验与临床研究txt、pdf、电子书免费下载,MMP-3、MMP-9在创伤性深静脉血栓形成中作用的实验与临床研究在线阅读大结局--361学术--提供学术界所需的全部图书
法：本研究分为三个部分：第一部分：建立创伤性DVT大鼠模型,基因芯片检测大鼠股静脉组织基因表达变化,采用倍数变化分析法(Fold Change,FC)筛选出差异表达基因2458个其中上调1146个；下调1312个。设定参数Single Log2 ra
基于1个网页-
生物芯片术语 - 豆丁网
通过电泳分离蛋白混合物，然后再印迹法转移蛋白到硝化纤维等支持物，用 于蛋白和多肽检测的实验方法。 倍性变化（fold change ，FC）一种用于描述两个用于相比的对象数量差异的方法。例如，第一个样本和第二个样本的量是50/10，那么FC
基于1个网页-
倍数变化分析
皱襞的变化
Fold change
Fold change is a measure describing how much a quantity changes going from an initial to a final value.
For example, an initial value of 30 and a final value of 60 corresponds to a fold change of 2, or in common terms, a two-fold increase.
以上来源于:
Compared with the placebo group, on average, participants who took paroxetine had a 6.8-fold change in signs of neuroticism, and 3.5 times as much change in signs of extroversion.
与服用安慰剂的小组相比，服用帕罗西汀的参与者在神经过敏症状的变化是前者的6.8倍，而在外向性的变化上则是前者的3.5倍。
Compared with the placebo group, on average, participants who took paroxetine had a 6.8-fold change in signs of neuroticism, and 3.5 times as much change in signs of extroversion.
与服用安慰剂的对照组相比，服用帕罗西汀的参与者在神经过敏症状的变化是前者的6.8倍，而在外向性的变化上则是前者的3.5倍。
Folds of rock layers that have similar mechanical properties or competence tend to be harmonic, with little change in fold shape, symmetry or wavelength from one layer to the next.
相似力学特性质或强岩层的褶皱趋于成为谐调的褶皱，在褶皱形态上、对称性上、或者从一层到下一层的褶皱波长上没有什么变化。
The results were so convincing that the programme was expanded 100-fold despite a change of government.
$firstVoiceSent
- 来自原声例句
请问您想要如何调整此模块？
感谢您的反馈，我们会尽快进行适当修改！
请问您想要如何调整此模块？
感谢您的反馈，我们会尽快进行适当修改！浙江大学干细胞与转化医学联盟
&&& 为进一步促进浙江省和浙江大学干细胞与再生医学领域的创新、交流与合作，浙江大学干细胞联盟应运而生。日，其在浙江大学紫金港校区召开第一次研讨会，并宣布正式成立。
&&& 为提供一个良好的交流平台，联盟将根据各成员的背景和特长，定期举行教授沙龙活动，届时将通过官方网站和官方微信公众号等方式通知大家。
当前位置： >
数据分析-【差异基因list做KEGG】
时间：日 访问次数:172
差异基因list做KEGG
& 为扩大浙江大学干细胞与转化医学联盟的影响力，小编Andyhigh亲测可用如下代码，可做芯片或高通量测序数据得到的差异基因分析，并进行了简单注释。可针对所有基因，上调基因，下调基因做KEGG的富集分析。代码参考生信菜鸟团博主相关代码。
链接：http://www.bio-/
setwd(&C:/Users/Administrator/Desktop/R/CLL&)
rm(list=ls())
library(CLL)
data(sCLLex)
library(limma)
exprSet = exprs(sCLLex)
pdata=pData(sCLLex)
group_list = pdata$Disease
design=model.matrix(~factor(group_list))
fit=lmFit(exprSet,design)
fit=eBayes(fit)
topTable(fit,coef=2,n=20)
DEG&-topTable(fit,coef=2,n=20)&& ###此外20为选择前20个基因，如需要所有可改为&n=Inf&
DEG$symbol &- rownames(DEG)&&& #将DEG中的symbol归结为左侧栏目
write.table(DEG, file=&CLL.xls&, sep=&\t&,quote=F)
probeset=rownames(DEG[abs(DEG[,1])&0.5 & DEG[,4]&0.05,])
###上述为获取差异表达基因list；筛选标准为logFoldchange&0.5(包括上调和下调)，P&0.05
library(annotate)
library(hgu95av2.db)
platformDB=&hgu95av2.db&
EGID &- as.numeric(lookUp(probeset, platformDB, &ENTREZID&))
length(unique(EGID))
###将探针ID转变为EG ID，切EGID在第一列；
diff_gene_list &- unique(EGID)
diff_gene_list
source(&http://bioconductor.org/biocLite.R&)
options(BioC_mirror=&http://mirrors./bioc/&)
biocLite(&GOstats&)
biocLite(&org.Hs.eg.db&)
biocLite(&KEGG.db&)
library(GOstats)
library(org.Hs.eg.db)
library(KEGG.db)
#then do kegg pathway enrichment
hyperG.params = new(&KEGGHyperGParams&, geneIds=diff_gene_list, universeGeneIds=NULL, annotation=&org.Hs.eg.db&,
&&&&&&&&&&&&&&&&&&&&
&&&&&&&&&&&&&&&&&&&& categoryName=&KEGG&, pvalueCutoff=1, testDirection = &over&)
KEGG.hyperG.results = hyperGTest(hyperG.params)
htmlReport(KEGG.hyperG.results, file=&kegg.enrichment.html&, summary.args=list(&htmlLinks&=TRUE))
当然也可以用自己的数据进行分析，下面是相关代码
rm(list=ls())
setwd(&F:/R/2D3D/Log2/LX&)
eset&-read.table(&lx 2D3D.csv&, header=TRUE, row.names=&GeneID&, sep = &,&, )
group_list&-read.csv(&2D3D group.csv&, row.names=&SampleID&, sep=&,&)
group &- group_list$System
library(limma)
design=model.matrix(~factor(group))&& &
fit=lmFit(eset,design)&& #将表达量和design 线性拟合
fit=eBayes(fit)&&&&&&&&&&&&&&&&&& #利用Bayes方法进行检验
options(digits = 2)
DEG &- topTable(fit,coef=2,adjust='BH',n=Inf)
DEG$symbol &- rownames(DEG)
write.table(DEG, file=&lx 2D3D.xls&, sep=&\t&,quote=F)
source(&http://bioconductor.org/biocLite.R&)
options(BioC_mirror=&http://mirrors./bioc/&)
biocLite(&GOstats&)
biocLite(&org.Hs.eg.db&)
biocLite(&KEGG.db&)
library(GOstats)
library(org.Hs.eg.db)
library(KEGG.db)
probeset=rownames(DEG[abs(DEG[,1])&0.5 & DEG[,4]&0.05,])
length(unique(probeset))
diff_gene_list &- unique(probeset)
diff_gene_list
hyperG.params = new(&KEGGHyperGParams&, geneIds=diff_gene_list, universeGeneIds=NULL, annotation=&org.Hs.eg.db&,
&&&&&&&&&&&&&&&&&& &
&&&&&&&&&&&&&&&&&&& categoryName=&KEGG&, pvalueCutoff=1, testDirection = &over&)
KEGG.hyperG.results = hyperGTest(hyperG.params)
htmlReport(KEGG.hyperG.results, file=&kegg.enrichment2.html&, summary.args=list(&htmlLinks&=TRUE))