Skip to content

CNVkit for RNA-seq: log2 in *.cnr file is all zero #847

Description

@jingyu9603

I am trying to analyze RNA expression data using CNVkit. Here is the code I used:
cnvkit.py import-rna ./SRR8613811.txt --gene-resource ../ensembl-gene-info.hg38.tsv --correlations ~/softwave/cnvkit-0.9.10/data/tcga-skcm.cnv-expr-corr.tsv --output test2-summary.tsv --output-dir test2/
The sample read count file has the following form:

Ensembl_ID SRR8613811
ENSG00000000003 365
ENSG00000000005 8
ENSG00000000419 691
ENSG00000000457 711
ENSG00000000460 602
ENSG00000000938 94
ENSG00000000971 820
ENSG00000001036 566
ENSG00000001084 719
ENSG00000001167 1575
ENSG00000001460 320
ENSG00000001461 2061

The code runs smoothly

Dropping 162 / 25299 rarely expressed genes from input samples
Loading gene metadata and TCGA gene expression/CNV profiles
Loaded ../ensembl-gene-info.hg38.tsv with shape: (221323, 9)
│Loaded /home/yj2/softwave/cnvkit-0.9.10/data/tcga-skcm.cnv-expr-corr.tsv with shape: (19177, 4)
Resetting 2846 ambiguous genes' correlation coefficients to default 0.100000
Trimmed gene info table to shape: (63966, 13)
Aligning gene info to sample gene counts
Weighting genes with below-average read counts
Calculating normalized gene read depths
Weighting genes by spread of read depths
Weighting genes by spearman_r correlation coefficient
Weighting genes by pearson_r correlation coefficient
Weighting genes by kendall_t correlation coefficient
/home/yj2/miniconda3/lib/python3.7/site-packages/scipy/stats/stats.py:283: RuntimeWarning: divide by zero encountered in log
log_a = np.log(a)
Writing output files
Wrote test2-summary.tsv with 24951 rows
Wrote test2/SRR8613811.cnr with 24951 regions

However, when I checked the *.cnr file, I found that the log2 of all the genes was 0. AS follows:

chromosome start end gene log2 depth gc tx_length weight
1 11869 14409 DDX11L1 0 0.54315 0.575 1657 1
1 14404 29570 WASH7P 0 9.10437 0.5443 1351 1
1 131025 134836 CICP27 0 6.29591 0.5971 3812 1
1 185217 195411 FO538757.1 0 10.952 0.5654 1397 1
1 257864 359681 AP006222.1 0 23.6302 0.3998 5603 1
1 347982 348366 RPL23AP24 0 0.519481 0.5039 385 1
1 365389 522928 AL732372.2 0 148.559 0.4261 902 1
1 450703 451697 OR4F29 0 0.301508 0.4573 995 1
1 491225 493241 AL732372.3 0 2.82486 0.6619 1239 1

I don't know what went wrong. Can someone help me? Thank you very much!

Metadata

Metadata

Assignees

No one assigned

    Labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions