搜索
查看: 2623|回复: 3

[其他] TCGA分析遇到的报错,请帮忙看看,谢谢

[复制链接]

4

主题

20

帖子

395

积分

中级会员

Rank: 3Rank: 3

积分
395
发表于 2017-3-28 22:16:13 | 显示全部楼层 |阅读模式
尊敬的群主和大神们,
我在练习“三分钟利用TCGA数据发表文章中的图表”,(链接为:http://3g.sanwen.net/a/ddatkpo.html)时遇到了一个问题,前面的命令都能顺利执行
但到了命令“
while(all.Found == F){
  mRNA.Exp[[page.Counter]] = Samples.mRNASeq(format = "csv",
                                             gene = diff.Exp.Genes,
                                             cohort = "PAAD",
                                             tcga_participant_barcode =
                                             paad.Pats$tcga_participant_barcode,
                                             page_size = page.Size,
                                             page = page.Counter)
  if(nrow(mRNA.Exp[[page.Counter]]) < page.Size)
    all.Found = T
  else
    page.Counter = page.Counter + 1
}

这步时,报错为“Error in mRNA.Exp[[page.Counter]] : subscript out of bounds
In addition: Warning message:
In download.Data(url, format, page) :
  The API responded with code  504. Your query might be to big”
百度和谷歌也没有找到解决方法,我也通过该帖子的微信号请教,但仍没有回复,
所以想请教群主和大神们帮忙解决,谢谢!

============================

以下是该帖子的全部代码:
install.packages("devtools")
install.packages("ggplot2")
devtools::install_github("mariodeng/FirebrowseR")
require(FirebrowseR)

cohorts = Metadata.Cohorts(format = "csv")
cancer.Type = cohorts[grep("Pancreatic adenocarcinoma", cohorts$description, ignore.case = T), 1]
print(cancer.Type)

all.Received = F
page.Counter = 1
page.size = 150
paad.Pats = list()

while(all.Received == F){
  paad.Pats[[page.Counter]] = Samples.Clinical(format = "csv",
                                               cohort = cancer.Type,
                                               page_size = page.size,
                                               page = page.Counter)

  if(page.Counter > 1)
    colnames(paad.Pats[[page.Counter]]) = colnames(paad.Pats[[page.Counter-1]])
  if(nrow(paad.Pats[[page.Counter]]) < page.size){
    all.Received = T

  } else{
    page.Counter = page.Counter + 1
  }}

paad.Pats = do.call(rbind, paad.Pats)
dim(paad.Pats)


diff.Exp.Genes = c("ESR1", "GATA3", "XBP1", "FOXA1", "ERBB2", "GRB7", "EGFR","FOXC1", "MYC")
all.Found = F
page.Counter = 1
mRNA.Exp = list()
page.Size = 2000 # using a bigger page size is faster
while(all.Found == F){
  mRNA.Exp[[page.Counter]] = Samples.mRNASeq(format = "csv",
                                             gene = diff.Exp.Genes,
                                             cohort = "PAAD",
                                             tcga_participant_barcode =                                             

                                               paad.Pats$tcga_participant_barcode,
                                             page_size = page.Size,
                                             page = page.Counter)
  if(nrow(mRNA.Exp[[page.Counter]]) < page.Size)
    all.Found = T
  else
    page.Counter = page.Counter + 1
}
mRNA.Exp = do.call(rbind, mRNA.Exp)
dim(mRNA.Exp)
mRNA.Exp = mRNA.Exp[which(mRNA.Exp$sample_type %in% c("NT", "TP")), ]
mRNA.Exp$expression_log2=as.numeric(mRNA.Exp$expression_log2)


library(ggplot2)

p = ggplot(mRNA.Exp, aes(factor(gene), z.score))
p +
  geom_boxplot(aes(fill = factor(sample_type))) +
  scale_y_continuous(limits = c(-2, 14)) +
  scale_fill_discrete(name = "Tissue")


library(ggplot2)
p = ggplot(mRNA.Exp, aes(factor(gene), expression_log2))
p +
  geom_boxplot(aes(fill = factor(sample_type))) +
  scale_y_continuous(limits = c(0, 8)) +
  scale_fill_discrete(name = "Tissue")



write.table(mRNA.Exp,file="C:/Users/Administrator/Desktop/mRNA.xls",sep="\t",quote=F)



回复

使用道具 举报

4

主题

24

帖子

184

积分

注册会员

Rank: 2

积分
184
发表于 2017-3-29 03:47:06 | 显示全部楼层
本帖最后由 tsznxx 于 2017-3-29 03:52 编辑

504错误是网页响应时间过久。应该是patient 太多了,时间响应太久。把paad.Pats$tcga_participant_barcode改成paad.Pats$tcga_participant_barcode[1:10]就可以通过。但是你要下载所有病人的话,就比较麻烦了。

这个是FireBrowseR里下载用的URL,你看着改改,可以直接得到csv文件。
[AppleScript] 纯文本查看 复制代码
http://firebrowse.org/api/v1/Samples/mRNASeq?format=csv&gene=ESR1%2CGATA3%2CXBP1%2CFOXA1%2CERBB2%2CGRB7%2CEGFR%2CFOXC1%2CMYC&cohort=PAAD&tcga_participant_barcode=TCGA-IB-7652%2CTCGA-FB-A4P6%2CTCGA-HZ-7289%2CTCGA-FB-AAQ6%2CTCGA-IB-7646&protocol=RSEM&page=1&page_size=2000&sort_by=cohort


PS: 这些东西感觉华而不实,看似捷径,花狸狐哨的整了一堆,其实都是黑盒子,不知道里面到底在弄啥,耽误很多时间,也学不到什么东西。就算跑通了,哪天R/firebrowseR/FireBrowse网站一更新,又不知道会出什么鬼问题。
Firehose有下载链接,直接下载全部基因的表达量不就完了,数据是表格形式的:geneXsampleID。关于sampleID,可以搜索下 TCGA Analyte ID,里面有解释。然后就是把你要的基因挑出来,按照cancertype分组,做个图罢了。(注:URL里改改cancertype,其他的一并都下载了)
[AppleScript] 纯文本查看 复制代码
http://gdac.broadinstitute.org/runs/stddata__2016_01_28/data/PAAD/20160128/gdac.broadinstitute.org_PAAD.Merge_rnaseqv2__illuminahiseq_rnaseqv2__unc_edu__Level_3__RSEM_genes_normalized__data.Level_3.2016012800.0.0.tar.gz


如果非要某个基因的表达量,可以用Python 的highlevel API(支持命令行调用)下载数据,一行命令搞定。
https://confluence.broadinstitute.org/display/GDAC/fbget

[Python] 纯文本查看 复制代码
import fbget
print fbget.mrnaseq("egfr", cohort="ucs")
tcga_participant_barcode        gene        expression_log2        z-score        cohort        sample_type        protocol        geneID
TCGA-QN-A5NN        EGFR        7.06162500905        -0.59899352506        UCS        TP        RSEM        1956
TCGA-QM-A5NM        EGFR        8.16734387649        -0.298443593752        UCS        TP        RSEM        1956
TCGA-NG-A4VW        EGFR        8.93092623547        0.0932667888031        UCS        TP        RSEM        1956


至少你清楚的看到你得到的是什么东西。后面接着往下做就是了。




回复 支持 5 反对 0

使用道具 举报

3

主题

43

帖子

212

积分

中级会员

Rank: 3Rank: 3

积分
212
发表于 2017-3-29 08:52:43 | 显示全部楼层
tsznxx 发表于 2017-3-29 03:47
504错误是网页响应时间过久。应该是patient 太多了,时间响应太久。把paad.Pats$tcga_participant_barcode ...

赞一个~~~~
回复 支持 反对

使用道具 举报

4

主题

20

帖子

395

积分

中级会员

Rank: 3Rank: 3

积分
395
 楼主| 发表于 2017-3-29 09:07:07 | 显示全部楼层
谢谢“tsznxx”的回复,这个回复对我很有启发!赞!
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-10-22 23:01 , Processed in 0.033375 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.