Dear BioCommunity,

i recently implemented the R package TCGAbiolinks, to download raw HTSEQ counts for a provinsional cancer TCGA dataset (COAD). My code used was the following:

query.exp.hg38 <- GDCquery(project = "TCGA-COAD", 
data.category = "Transcriptome Profiling", 
data.type = "Gene Expression Quantification", 
workflow.type = "HTSeq - Counts")

GDCdownload(query.exp.hg38,files.per.chunk = 50)

exp.hg38 <- GDCprepare(query = query.exp.hg38,
save=TRUE,directory=mypath,save.filename = "exp.coad.htseqcounts.rda")

exp.hg38
class: RangedSummarizedExperiment 
dim: 56963 521 
metadata(1): data_release
assays(1): HTSeq - Counts
rownames(56963): ENSG00000000003 ENSG00000000005 ...
  ENSG00000281912 ENSG00000281920
rowData names(3): ensembl_gene_id external_gene_name
  original_ensembl_gene_id
colnames(521): TCGA-3L-AA1B-01A-11R-A37K-07
  TCGA-DM-A1D8-01A-11R-A155-07 ...
  TCGA-AA-3675-01A-02R-0905-07 TCGA-G4-6323-01A-11R-1723-07
colData names(101): sample patient ...
  subtype_vascular_invasion_present subtype_vital_status

test <- TCGAanalyze_Normalization(tabDF=exp.hg38,geneInfo)
I Need about  0 seconds for this Complete Normalization Upper Quantile  [Processing 80k elements /s]  
Step 1 of 4: newSeqExpressionSet ...
Step 2 of 4: withinLaneNormalization ...
Error in names(y) <- 1:length(y) : 
  'names' attribute [2] must be the same length as the vector [0]
Timing stopped at: 0 0 0.04 

 test2 <- TCGAbiolinks::TCGAanalyze_Normalization(tabDF = exp.hg38,geneInfo=TCGAbiolinks::geneInfoHT,method="gcContent")
I Need about  371 seconds for this Complete Normalization Upper Quantile  [Processing 80k elements /s]  
Step 1 of 4: newSeqExpressionSet ...
Step 2 of 4: withinLaneNormalization ...
Step 3 of 4: betweenLaneNormalization ...
Step 4 of 4: .quantileNormalization ...
> dim(test3)
[1] 23238   521

Any ideas about the first above error ? as also why the second approach worked,  but the different number  of features in the second approach is reduced to more than the half of the 56963 initial number ? From a small search in the function, in mentions in the geneInfo argument:

geneInfo= Information matrix of 20531 genes about geneLength and gcContent. Two objects are provided: TCGAbiolinks::geneInfoHT,TCGAbiolinks::geneInfo

But, again how it explains the above issue ? perhaps it has to do something with the EDASeq R package, which is mentioned for implementing the above normalization function ?

Briefly, my main goals here, is to perform downstream DE analysis between cancer and normal samples, as also a suarvival analysis with a specific subset of genes.

Thank you in advance,

Efstathios

Dear @svlachavas (Efstathios) thank you for interest of using TCGAbiolinks for your research. I think that your post it is similar to the one I answered 1 week ago in our TCGAbiolinks github issues here: https://github.com/BioinformaticsFMRP/TCGAbiolinks/issues/164. I am sorry that I wasn't reading also the bioconductor.support. Was my answer there useful? Thank you davide risso  (Davide) for helping us to find a solution. Again if you are using data from Gene expression aligned against hg38 in that case you can use geneInfoHT object in TCGAanalyze_Normalization as well documented in our vignette 

http://bioconductor.org/packages/release/bioc/vignettes/TCGAbiolinks/inst/doc/tcgaBiolinks.html#htseq_data:_downstream_analysis_brca

And thank you for your suggestion, I used the function getGeneLengthAndGCContent from Davide's package EDASeq, and when it possibile I will add a function to update the geneInfoHT in TCGAbiolinks.

Hi,

I am not familiar with the implementation of TCGAbiolinks and with how they use the functions in the EDASeq package, so I'm not fully able to address your problem. Hopefully, some of the TCGAbiolinks authors / maintainers will chip in.

The error seems related to some naming issues. How does your geneInfo object look like? Can you please post here the results of

class(geneInfo)
head(geneInfo)

I'm just guessing, but could it be that the geneInfo object uses a different set of IDs for genes compared to the exp.hg38 object? (e.g. Ensembl vs gene symbols?)

Also, from the help page that you are providing above, it looks like the TCGAbiolinks package only provides the gene length and GC content for 20531 genes. I'm not sure why you end up with 23238 genes. This is a question for the TCGAbiolinks authors.

Best, Davide