TCGA数据下载方式小总结

之前对TCGA做了简单的了解,粗略了解了什么是TCGA,TCGA是做什么的等情况,接下来肯定是要学会如何下载TCGA数据,毕竟只有下载了数据才能接下去继续实战学习

官网常规下载

TCGA自2016年改版后,下载方式变得大为不同,数据都整合在GDC(Genomic Data Commons)的DATA PORTAL中,所以我们要先进网址https://portal.gdc.cancer.gov/

我以乳腺癌(BRCA)为例,进入TCGA-BRCA界面,这次选择下载的是RNA-Seq数据,点击RNA-Seq的1092个case,然后在点击Files自行设定筛选条件(我选择下载count数据,记得还要选open的数据),最后加所有样本都加入购物车Cart(总共会有1222个文件)。为什么会比case数多呢,说明一些case不止对应1个RNA-Seq样本的count文件,有些是重复。

接下来就稍微有点不同了,如果只是比较少量的样本,直接点击Download即可下载,每个样本都是一个压缩文件。但是,TCGA数据库在数据下载有规定:让Cart文件夹大于50M时(这个依据网络情况,和下载用户数目),只能通过Data Transfer Tool工具进行下载。所以我这次要使用Data Transfer Tool工具来下载数据。

首先是要安装Data Transfer Tool,总共工具很简单,我是选择windows系统来安装的,下载后解压缩即可使用。至于怎么用,我刚开始看到是.exe文件,就以为是软件界面的,后来才发现原来是通过命令行下载的。先从刚才Cart界面点击Download选择下载Manifest,然后在windows的cmd中输入下载命令

gdc-client download -m gdc_manifest_20171102_125645.txt

然后就是慢慢下载了,如果中间有报错断了,一般就是网络不好的原因,重新下载或者换个时间段就行了,由于每一个下载的结果是一个样本的压缩文件,1222个样本就是1222个压缩包,所以我要解压缩这些文件并合并成1个矩阵文件,主要可以分为3步:

  1. 将样本名称改为TCGA-C8-A12Q-01这种形式,这种ID可以从Cart界面下载metadata文件中获取,在metadata文件中提取下载文件名和样本名的对应关系,最开始是用perl的,但是感觉perl会比较繁琐,所以改用R来处理。选择的是rjson这个包,也试了下jsonlite包,基本功能是差不多的

    library(rjson)
    
    data <- fromJSON(file = "metadata.cart.2017-10-30T08-39-11.222982.json")
    
    data2 <- sapply(data, function(x){
      file_name <- x$file_name
      file_id <- x$file_id
      submitter_id <- x$case[[1]]$samples[[1]]$submitter_id
      tmp <- c(file_name, file_id, submitter_id)
    })
    
    data2 <- t(data2)
    colnames(data2) <- c("filename", "fileid", "submitterid")
    
    write.table(data2, file = "sample2id.txt", sep = "\t", col.names = F, quote = F, row.names = F)
    

    输出的sample2id.txt文件中总共有3列,第一列为我们需要解压的文件名,第二列为我们下载下来的每个压缩包所在的文件夹的名称,第3列则是submitter id,就是我们需要的样品名

  2. 接下来要对下载的文件进行批量解压缩,并输出到file文件夹中,该shell脚本输入参数为上述的sampleid.txt文件。这个代码有个问题就是如果遇到样品名相同的两个样品的,后者文件会覆盖前者文件。当然只要将每个样本的count文件输出在各自文件夹中就不会有这个错误了(考虑到这样做的话,后面还是需要去掉重复样本,所以暂时不这样做了)

    #!/bin/bash
    
    cat $1 |while read line
    do
      arr=($line)
      filename=${arr[0]}
      folder=${arr[1]}
      submitterid=${arr[2]}
      gunzip -c ./${folder}/${filename} > ./file/${submitterid}.count
    done
    

    运行上述bash extract_files.sh ./sample2id.txt脚本后,1222个TCGA乳腺癌的样本变成了1217个样本(说明有5个样本的名称是重复的,所以这些文件被覆盖了),这是在file文件夹中有1217个样本的count文件

  3. 最后将每个样本的count数进行合并,这次我用perl来整合成一个count矩阵文件,并且只保留01A和11A结尾的样本(01A表示肿瘤细胞,11表示癌旁样本)。如果以B结尾的样品名,则表明是重复,其实还有C结尾的,这次都先选择去除

    #!/bin/perl -w
    use strict;
    
    my $path = shift @ARGV;
    
    opendir DIR, $path or die;
    my @dir = readdir DIR;
    
    my $header;
    my @sample;
    my %hash;
    my %rm_rep;
    
    foreach my $file (@dir) {
        if ($file =~ /^(\w+.*-\d+)[A-Z]\.count/) {
            my $sample_id = $1;
            next unless ($sample_id =~ /01$/ || $sample_id =~ /11$/);
            if ($rm_rep{$sample_id}) {
                next;
            }else{
                $rm_rep{$sample_id} = 1;
            }
            push @sample, $sample_id;
            $header .= "\t$sample_id";
    
            open my $fh, $path.$file or die;
            while (<$fh>) {
                chomp;
                next if ($_ =~ /^_/);
                my @array = split /\t/, $_;
                $hash{$array[0]} -> {$sample_id} = $array[1];
            }
            close $fh;
        }
    }
    
    print "$header\n";
    map{
        my $gene = $_;
        print "$gene";
        foreach my $file (@sample) {
            print "\t".$hash{$gene} -> {$file};
        }
        print "\n";
    }keys %hash;
    

    运行perl merged_count.pl ./file/,代码看起来有点繁琐,但是还是可以用的~~最后输出一个整合后的count矩阵文件,后续分析就是根据这个count数据来了

Firehose下载

第一种方式个人觉得比较可控,下载的数据的准确性也是有保证的。但是如果觉得麻烦的话,可以使用Broad Institute’s GDAC Firehose,该网站不仅提供TCGA数据的下载还有其分析后的数据的展示。这里我只用过其下载的功能,所以其他的自行去网站http://gdac.broadinstitute.org/逛逛吧

进入网站后,点击上方的Dashboards下的Standard Data,然后再点击网页上的乳腺癌(BRCA)对应的open按钮,结果会出现一系列文件可供下载,为什么更好的跟第一种方法进行比较,我选择了下载RNA-Seq的count数据对应的文件,可以选择下载gdac.broadinstitute.org_BRCA.Merge_rnaseq__illuminahiseq_rnaseq__unc_edu__Level_3__gene_expression__data.Level_3.2016012800.0.0.tar.gz,其实还有个rnaseqv2对应文件,刚开始我也是下载错了,后来查了下才发现前者和后者是有区别的,官方给的解释是:rnaseqv2相比rnaseq采用不同的算法(我看了下似乎是用RSEM进行定量的?因为rnaseq是采用htseq来定量的),但是呢rnaseqv2相比rnaseq样本数多一点。。。

RNASeq Version 2 is similar to RNASeq in that it uses sequencing data to determine gene expression levels. RNASeq Version 2 uses a different set of algorithms to determine the expression levels are the results are presented in a slightly different set of files.

下载后将文件夹解压,里面有不少文件,我先只拿count文件来看看,这里就显示出Firehose下载的好处了,不用像第一种方法那样将一个个样本合并,其下载结果就是合并好的数据,而且已经将Ensembl id 转化成gene symbol,很是方便

但是!我发现一个问题,就是Firehose下载下来的文件的样本数跟第一种方法下载后并不一致,粗略看了下,大概少了几百个样本。。就算是rnaseqv2也少了几个样本。然后我将两个下载方法获得的样本名都进行去重,然后比较发现有几个样本虽然在第二个方法中出现,但是在TCGA官网上并没有RNA-Seq的数据了!!!!这时我才想起来Firehose数据并不是跟TCGA官方数据更新同步的,现在下载的Firehose的数据是2016年初的。。。

所以Firehose下载操作确实方便,但是数据并不是时时更新开放下载,如果不在意样本数目的话,这个方式确实是一个比较方便的下载TCGA数据的方法

生信人开发的小工具

我也是听别人介绍了才知道生信人公众号做了一个下载TCGA数据的小工具,初步尝试了下,感觉不错,操作蛮简单的,并且不仅提供了数据下载功能,还有合并样本、临床信息还有ID转化的功能

具体操作还是看生信人的宣传文章把https://www.shengxin.ren/article/27

这个优点我觉得在于下载数据的实时性有保证,并且它也是从官网下载的,最主要它支持断点下载。。。就算你网再不稳定。。睡一觉总能下好了,比第一种下载方式方便(因为第一种网不好容易断,导致一些文件无法被正常下载),而且最后整合的样本数目(去重后)跟第一种方法一致,这让我心里也有了底。。。。

还有一个,该工具在使用ID转化功能后,能将lncRNA区分开。。。这个就满足好多人的需求了。。但是我也准备自己来写个代码做这件事。。比较下结果。。

本文出自于http://www.bioinfo-scrounger.com转载请注明出处