search

ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
105 stars 30 forks source link

下载TCGA所有癌症的maf文件做signature分析 #13

Closed ixxmu closed 4 years ago

ixxmu commented 4 years ago

https://mp.weixin.qq.com/s?__biz=MzAxMDkxODM1Ng==&mid=2247484960&idx=1&sn=80641a93f36800d791eca95f41d4f2f9&chksm=9b48469bac3fcf8d7736fe5e04794db90dbc6d13b3261efa2ba14bd15c78aa6433c5059fa4ea&scene=21#wechat_redirect

github-actions[bot] commented 4 years ago

下载TCGA所有癌症的maf文件做signature分析 by 生信技能树

才sanger研究所已经做好了这个分析，但是值得我们重复一下，效果如下：


TCGA所有癌症的mutation signature

首先TCGA所有癌症的maf文件

maf格式的mutation记录文件在TCGA里面已经是level4的数据啦，所以是完全open的，可以随意下载，只需要去其GDC官网简单点击，选择即可。

主要步骤就是在https://portal.gdc.cancer.gov/repository里面点击过滤文件类型，选择maf格式，再过滤access权限，选择open即可，最后得到的132个文件就是我们需要的。

总共是2.19GB的文件，每个癌症种类都有4种maf文件，分别是用mutect,muse,vanscan,somaticsniper这4款软件call 到的somatic mutation文件。

下载方式这里我选择下载它们132个文件的manifest文件，然后用GDC提供的官方工具来下载！关于这个工具，我 在生信技能树论坛写过教程，就不多说了，自己去看哈，现在下载TCGA数据也是非常方便，首先是GDC网站及客户端 就是安装成功后，运行 ./gdc-client download -m manifest_xxx.txt j即可。这个manifest文件就是自己刚才创造并且下载的。

  1. cd ~/institute/TCGA/GDC_NCBI/all 
  2. ~/biosoft/GDC/gdc-client download -m gdc_manifest.2017-08-25T02-57-11.281090.txt

但是这个工具，提供的电脑操作系统版本有限哦If you are a user of CentOS 6 or RedHat Enterprise Release 6 and wish to use the Data Transfer Tool, contact the GDC Help Desk for assistance.

所以我是在MAC里面下载好了，再上传到我的服务器去的！

然后根据MAF文件制作signature

我是根据这篇文章Mutational Profile of Metastatic Breast Cancers: A Retrospective Analysis里面的方法来做的，他们的方法描述如下：

  1. De novo mutational signature analysis was done using the Matlab Welcome Trust Sanger Institute’s signature framework. We used the deconstructSigs R package to determine the contribution of the known signatures that explain each sample mutational profile with more than 50 somatic mutations. We considered the 13 signatures (Signatures 1, 2, 3, 5, 6, 8, 10, 13, 17, 18, 20, 26, and 30) operative in breast cancer as defined in COSMIC (http://cancer.sanger.ac.uk/signatures/matrix.png). A signature was defined as operative or predominant if its contribution to the mutational pattern was respectively >25% (or >100 mutations) or >50%.

虽然我以前写过一个类似的教程，用SomaticSignatures包来解析maf突变数据获得mutation signature 这里还是再学习学习这个新的工具deconstructSigs R package吧。

这是个R包，所以直接在Rstudio里面安装即可，这里选取BRCA的somatic mutation的MAF文件做一下分析，看看四个软件找出的变异，是否在signature上面有差异。

  1. install.packages('deconstructSigs')
  2. # dependencies 'BSgenome', 'BSgenome.Hsapiens.UCSC.hg19' 
  3. BiocInstaller::biocLite('BSgenome')
  4. BiocInstaller::biocLite('BSgenome.Hsapiens.UCSC.hg19')
  5. ## https://github.com/raerose01/deconstructSigs
  7. file1='TCGA.STAD.muse..somatic.maf.gz'
  8. TCGA.STAD.muse=read.table(file1,sep = '\t',quote="",header = T)
  9. TCGA.STAD.muse[1:5,1:15]
  10. ## data frame including 5 columns: sample.ID,chr,pos,ref,alt 
  11. sample.mut.ref <- data.fram(Sample='TCGA.STAD.muse', 
  12.                                 chr = TCGA.STAD.muse[,5], 
  13.                                 pos = TCGA.STAD.muse[,6], 
  14.                                 ref = TCGA.STAD.muse[,11], 
  15.                                 alt = TCGA.STAD.muse[,13])
  16. sigs.input <- mut.to.sigs.input(mut.ref = sample.mut.ref, 
  17.                                 sample.id = "Sample", 
  18.                                 chr = "chr", 
  19.                                 pos = "pos", 
  20.                                 ref = "ref", 
  21.                                 alt = "alt")
  22. class(sigs.input)
  24. sample_1 = whichSignatures(tumor.ref = sigs.input, 
  25.                            signatures.ref = signatures.nature2013, 
  26.                            sample.id = 'TCGA.STAD.muse', 
  27.                            contexts.needed = TRUE,
  28.                            tri.counts.method = 'exome')
  30. # Plot example                  
  31. plot_example <- whichSignatures(tumor.ref = sigs.input  ,  
  32.                        signatures.ref = signatures.nature2013, 
  33.                        sample.id = 'TCGA.STAD.muse' )
  35. # Plot output
  36. plotSignatures(plot_example, sub = 'example')

这个里面有一个问题，就是deconstructSigs R package似乎只支持hg19版本的基因组，而我下载的TCGA的MAF是hg38版本的，所以代码虽然是对的，但实际上做出的结果是不对的，需要把下载的TCGA的maf文件进行坐标转换。

而所谓批量，无非就是在上面的R脚本里面加入一个循环咯。

（点击阅读原文有这个包的详细说明书哈！）

注意事项，下载的MAF文件可能有两种格式 ，可能是47列，或者120列，第一行一般都是 头文件，注释着每一列的信息，的确，信息量有点略大。如下：

  1.      1    Hugo_Symbol
  2.      2    Entrez_Gene_Id
  3.      3    Center
  4.      4    NCBI_Build
  5.      5    Chromosome
  6.      6    Start_Position
  7.      7    End_Position
  8.      8    Strand
  9.      9    Consequence
  10.     10    Variant_Classification
  11.     11    Variant_Type
  12.     12    Reference_Allele
  13.     13    Tumor_Seq_Allele1
  14.     14    Tumor_Seq_Allele2
  15.     15    dbSNP_RS
  16.     16    dbSNP_Val_Status
  17.     17    Tumor_Sample_Barcode
  18.     18    Matched_Norm_Sample_Barcode
  19.     19    Match_Norm_Seq_Allele1
  20.     20    Match_Norm_Seq_Allele2
  21.     21    Tumor_Validation_Allele1
  22.     22    Tumor_Validation_Allele2
  23.     23    Match_Norm_Validation_Allele1
  24.     24    Match_Norm_Validation_Allele2
  25.     25    Verification_Status
  26.     26    Validation_Status
  27.     27    Mutation_Status
  28.     28    Sequencing_Phase
  29.     29    Sequence_Source
  30.     30    Validation_Method
  31.     31    Score
  32.     32    BAM_File
  33.     33    Sequencer
  34.     34    t_ref_count
  35.     35    t_alt_count
  36.     36    n_ref_count
  37.     37    n_alt_count
  38.     38    HGVSc
  39.     39    HGVSp
  40.     40    HGVSp_Short
  41.     41    Transcript_ID
  42.     42    RefSeq
  43.     43    Protein_position
  44.     44    Codons
  45.     45    Hotspot
  46.     46    cDNA_change
  47.     47    Amino_Acid_Change
  1.      1    Hugo_Symbol
  2.      2    Entrez_Gene_Id
  3.      3    Center
  4.      4    NCBI_Build
  5.      5    Chromosome
  6.      6    Start_Position
  7.      7    End_Position
  8.      8    Strand
  9.      9    Variant_Classification
  10.     10    Variant_Type
  11.     11    Reference_Allele
  12.     12    Tumor_Seq_Allele1
  13.     13    Tumor_Seq_Allele2
  14.     14    dbSNP_RS
  15.     15    dbSNP_Val_Status
  16.     16    Tumor_Sample_Barcode
  17.     17    Matched_Norm_Sample_Barcode
  18.     18    Match_Norm_Seq_Allele1
  19.     19    Match_Norm_Seq_Allele2
  20.     20    Tumor_Validation_Allele1
  21.     21    Tumor_Validation_Allele2
  22.     22    Match_Norm_Validation_Allele1
  23.     23    Match_Norm_Validation_Allele2
  24.     24    Verification_Status
  25.     25    Validation_Status
  26.     26    Mutation_Status
  27.     27    Sequencing_Phase
  28.     28    Sequence_Source
  29.     29    Validation_Method
  30.     30    Score
  31.     31    BAM_File
  32.     32    Sequencer
  33.     33    Tumor_Sample_UUID
  34.     34    Matched_Norm_Sample_UUID
  35.     35    HGVSc
  36.     36    HGVSp
  37.     37    HGVSp_Short
  38.     38    Transcript_ID
  39.     39    Exon_Number
  40.     40    t_depth
  41.     41    t_ref_count
  42.     42    t_alt_count
  43.     43    n_depth
  44.     44    n_ref_count
  45.     45    n_alt_count
  46.     46    all_effects
  47.     47    Allele
  48.     48    Gene
  49.     49    Feature
  50.     50    Feature_type
  51.     51    One_Consequence
  52.     52    Consequence
  53.     53    cDNA_position
  54.     54    CDS_position
  55.     55    Protein_position
  56.     56    Amino_acids
  57.     57    Codons
  58.     58    Existing_variation
  59.     59    ALLELE_NUM
  60.     60    DISTANCE
  61.     61    TRANSCRIPT_STRAND
  62.     62    SYMBOL
  63.     63    SYMBOL_SOURCE
  64.     64    HGNC_ID
  65.     65    BIOTYPE
  66.     66    CANONICAL
  67.     67    CCDS
  68.     68    ENSP
  69.     69    SWISSPROT
  70.     70    TREMBL
  71.     71    UNIPARC
  72.     72    RefSeq
  73.     73    SIFT
  74.     74    PolyPhen
  75.     75    EXON
  76.     76    INTRON
  77.     77    DOMAINS
  78.     78    GMAF
  79.     79    AFR_MAF
  80.     80    AMR_MAF
  81.     81    ASN_MAF
  82.     82    EAS_MAF
  83.     83    EUR_MAF
  84.     84    SAS_MAF
  85.     85    AA_MAF
  86.     86    EA_MAF
  87.     87    CLIN_SIG
  88.     88    SOMATIC
  89.     89    PUBMED
  90.     90    MOTIF_NAME
  91.     91    MOTIF_POS
  92.     92    HIGH_INF_POS
  93.     93    MOTIF_SCORE_CHANGE
  94.     94    IMPACT
  95.     95    PICK
  96.     96    VARIANT_CLASS
  97.     97    TSL
  98.     98    HGVS_OFFSET
  99.     99    PHENO
  100.    100    MINIMISED
  101.    101    ExAC_AF
  102.    102    ExAC_AF_Adj
  103.    103    ExAC_AF_AFR
  104.    104    ExAC_AF_AMR
  105.    105    ExAC_AF_EAS
  106.    106    ExAC_AF_FIN
  107.    107    ExAC_AF_NFE
  108.    108    ExAC_AF_OTH
  109.    109    ExAC_AF_SAS
  110.    110    GENE_PHENO
  111.    111    FILTER
  112.    112    CONTEXT
  113.    113    src_vcf_id
  114.    114    tumor_bam_uuid
  115.    115    normal_bam_uuid
  116.    116    case_id
  117.    117    GDC_FILTER
  118.    118    COSMIC
  119.    119    MC3_Overlap
  120.    120    GDC_Validation_Status

猜你喜欢

基因组 | 游记 | 工作资讯 

学习课程 | 好书分享 


菜鸟入门

Linux | Perl | R语言 | 可视化 

R包 | perl模块 | python模块


数据分析

ChIP-seq（上）| ChIP-seq（下）| RNA-seq | miRNA

WGS,WES,RNA-seq组与ChIP-seq之间的异同


编程实践

第0题 | 探索人类基因组序列 | 最后一题


直播基因组分析

我的基因组 | 解惑帖

一个标准的基因检测报告目录

生信技能树

想真正看到R包用法，那么阅读原文是必须的，如果你比较拮据，就不需要赞赏啦，文末有广告，我猜你应该知道该怎么做的。