a-slide / pycoMeth

DNA methylation analysis downstream to Nanopolish for Oxford Nanopore DNA sequencing datasets
GNU General Public License v3.0
34 stars 2 forks source link

pycoMeth.common.pycoMethError: Invalid input file type passed (nanopolish_fn). Expecting Nanopolish call_methylation output TSV file #49

Closed jolespin closed 3 years ago

jolespin commented 3 years ago

I tried running the following command:

/usr/local/devel/ANNOTATION/jespinoz/anaconda3/envs/ont_env/bin/pycoMeth CpG_Aggregate --nanopolish_fn /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_3_date_8_3-20201019_P119_JIRA_1521_Batch7_sheared-barcode05/intermediate/call_output/methylation_calls.cpg.tsv.gz /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_6_date_8_3-20201021_P122_JIRA_1521_Batch9-barcode12/intermediate/call_output/methylation_calls.cpg.tsv.gz /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_2_date_8_3-20200916_P111_JIRA_1521_Batch2-barcode06/intermediate/call_output/methylation_calls.cpg.tsv.gz /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_1_date_8_3-20200908_P110_JIRA_1521_Batch1-barcode02/intermediate/call_output/methylation_calls.cpg.tsv.gz /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_6_date_8_3-20200916_P111_JIRA_1521_Batch2-barcode08/intermediate/call_output/methylation_calls.cpg.tsv.gz /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_1_date_8_3-20200918_P115_JIRA_1521_Batch3_reload-barcode09/intermediate/call_output/methylation_calls.cpg.tsv.gz /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_3_date_8_3-20200918_P115_JIRA_1521_Batch3_reload-barcode10/intermediate/call_output/methylation_calls.cpg.tsv.gz /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_1_date_8_3-20201207_P125_JIRA_1521_finalbatch_shear-barcode20/intermediate/call_output/methylation_calls.cpg.tsv.gz /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_1_date_8_3-20201023_P123_JIRA_1521_Batch10_sheared-barcode20/intermediate/call_output/methylation_calls.cpg.tsv.gz /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_5_date_8_3-20200928_P117_JIRA_1521_Batch5-barcode18/intermediate/call_output/methylation_calls.cpg.tsv.gz /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_2_date_8_3-20201008_P118_JIRA_1521_Batch6_sheared-barcode02/intermediate/call_output/methylation_calls.cpg.tsv.gz /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_5_date_8_3-20201207_P125_JIRA_1521_finalbatch_shear-barcode18/intermediate/call_output/methylation_calls.cpg.tsv.gz --ref_fasta_fn /usr/local/projdata/8570/projects/EPIENV3/jespinoz/Genomes/Rnor_6.0/Rattus_norvegicus.Rnor_6.0.dna_rm.toplevel.fa --output_bed_fn pycometh_output/treatment_8-3/intermediate/cpg_aggregate_output/CpG_aggregate.output.bed --output_tsv_fn pycometh_output/treatment_8-3/intermediate/cpg_aggregate_output/CpG_aggregate.output.tsv --sample_id treatment_8-3 --verbose --progress

But I got the following error. Is there a different syntax I need to use to include multiple files in the input?

Also, if I give a file that includes a list of file names, what should the header be?

Can these be gzipped? It says they can but I just want to double check.

(ont_env) -bash-4.1$ cat pycometh_output/treatment_8-3/log/1__cpg_aggregate.*
## Checking options and input files ##
    [DEBUG]: Options summary
    [DEBUG]:    Package name: pycoMeth
    [DEBUG]:    Package version: 0.4.25
    [DEBUG]:    Timestamp: 2021-01-22 20:31:13.061467
    [DEBUG]:    nanopolish_fn: ['/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_3_date_8_3-20201019_P119_JIRA_1521_Batch7_sheared-barcode05/intermediate/call_output/methylation_calls.cpg.tsv.gz', '/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_6_date_8_3-20201021_P122_JIRA_1521_Batch9-barcode12/intermediate/call_output/methylation_calls.cpg.tsv.gz', '/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_2_date_8_3-20200916_P111_JIRA_1521_Batch2-barcode06/intermediate/call_output/methylation_calls.cpg.tsv.gz', '/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_1_date_8_3-20200908_P110_JIRA_1521_Batch1-barcode02/intermediate/call_output/methylation_calls.cpg.tsv.gz', '/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_6_date_8_3-20200916_P111_JIRA_1521_Batch2-barcode08/intermediate/call_output/methylation_calls.cpg.tsv.gz', '/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_1_date_8_3-20200918_P115_JIRA_1521_Batch3_reload-barcode09/intermediate/call_output/methylation_calls.cpg.tsv.gz', '/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_3_date_8_3-20200918_P115_JIRA_1521_Batch3_reload-barcode10/intermediate/call_output/methylation_calls.cpg.tsv.gz', '/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_1_date_8_3-20201207_P125_JIRA_1521_finalbatch_shear-barcode20/intermediate/call_output/methylation_calls.cpg.tsv.gz', '/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_1_date_8_3-20201023_P123_JIRA_1521_Batch10_sheared-barcode20/intermediate/call_output/methylation_calls.cpg.tsv.gz', '/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_5_date_8_3-20200928_P117_JIRA_1521_Batch5-barcode18/intermediate/call_output/methylation_calls.cpg.tsv.gz', '/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_2_date_8_3-20201008_P118_JIRA_1521_Batch6_sheared-barcode02/intermediate/call_output/methylation_calls.cpg.tsv.gz', '/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_5_date_8_3-20201207_P125_JIRA_1521_finalbatch_shear-barcode18/intermediate/call_output/methylation_calls.cpg.tsv.gz']
    [DEBUG]:    ref_fasta_fn: /usr/local/projdata/8570/projects/EPIENV3/jespinoz/Genomes/Rnor_6.0/Rattus_norvegicus.Rnor_6.0.dna_rm.toplevel.fa
    [DEBUG]:    output_bed_fn: pycometh_output/treatment_8-3/intermediate/cpg_aggregate_output/CpG_aggregate.output.bed
    [DEBUG]:    output_tsv_fn: pycometh_output/treatment_8-3/intermediate/cpg_aggregate_output/CpG_aggregate.output.tsv
    [DEBUG]:    min_depth: 10
    [DEBUG]:    sample_id: treatment_8-3
    [DEBUG]:    min_llr: 2
    [DEBUG]:    verbose: True
    [DEBUG]:    quiet: False
    [DEBUG]:    progress: True
    [DEBUG]:    kwargs
    [DEBUG]:        subcommands: CpG_Aggregate
    [DEBUG]:        func: <function CpG_Aggregate at 0x2adbc168bf80>
## Parsing methylation_calls file ##
    [DEBUG]: Opening file /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_3_date_8_3-20201019_P119_JIRA_1521_Batch7_sheared-barcode05/intermediate/call_output/methylation_calls.cpg.tsv.gz in gzip mode
    [DEBUG]: Opening file /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_6_date_8_3-20201021_P122_JIRA_1521_Batch9-barcode12/intermediate/call_output/methylation_calls.cpg.tsv.gz in gzip mode
    [DEBUG]: Opening file /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_2_date_8_3-20200916_P111_JIRA_1521_Batch2-barcode06/intermediate/call_output/methylation_calls.cpg.tsv.gz in gzip mode
    [DEBUG]: Opening file /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_1_date_8_3-20200908_P110_JIRA_1521_Batch1-barcode02/intermediate/call_output/methylation_calls.cpg.tsv.gz in gzip mode
    [DEBUG]: Opening file /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_6_date_8_3-20200916_P111_JIRA_1521_Batch2-barcode08/intermediate/call_output/methylation_calls.cpg.tsv.gz in gzip mode
    [DEBUG]: Opening file /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_1_date_8_3-20200918_P115_JIRA_1521_Batch3_reload-barcode09/intermediate/call_output/methylation_calls.cpg.tsv.gz in gzip mode
    [DEBUG]: Opening file /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_3_date_8_3-20200918_P115_JIRA_1521_Batch3_reload-barcode10/intermediate/call_output/methylation_calls.cpg.tsv.gz in gzip mode
    [DEBUG]: Opening file /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_1_date_8_3-20201207_P125_JIRA_1521_finalbatch_shear-barcode20/intermediate/call_output/methylation_calls.cpg.tsv.gz in gzip mode
    [DEBUG]: Opening file /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_1_date_8_3-20201023_P123_JIRA_1521_Batch10_sheared-barcode20/intermediate/call_output/methylation_calls.cpg.tsv.gz in gzip mode
    [DEBUG]: Opening file /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_5_date_8_3-20200928_P117_JIRA_1521_Batch5-barcode18/intermediate/call_output/methylation_calls.cpg.tsv.gz in gzip mode
    [DEBUG]: Opening file /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_2_date_8_3-20201008_P118_JIRA_1521_Batch6_sheared-barcode02/intermediate/call_output/methylation_calls.cpg.tsv.gz in gzip mode
    [DEBUG]: Opening file /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_5_date_8_3-20201207_P125_JIRA_1521_finalbatch_shear-barcode18/intermediate/call_output/methylation_calls.cpg.tsv.gz in gzip mode
    [DEBUG]: Reading header from file: /local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_3_date_8_3-20201019_P119_JIRA_1521_Batch7_sheared-barcode05/intermediate/call_output/methylation_calls.cpg.tsv.gz
    [DEBUG]: Column names from header: 'chromosome / start / end / read_name / log_lik_ratio / log_lik_methylated / log_lik_unmethylated / num_calling_strands / num_cpgs / sequence'
    [DEBUG]: Input file type: unknown
    [DEBUG]: Closing file:/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_3_date_8_3-20201019_P119_JIRA_1521_Batch7_sheared-barcode05/intermediate/call_output/methylation_calls.cpg.tsv.gz
    [DEBUG]: Closing file:/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_6_date_8_3-20201021_P122_JIRA_1521_Batch9-barcode12/intermediate/call_output/methylation_calls.cpg.tsv.gz
    [DEBUG]: Closing file:/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_2_date_8_3-20200916_P111_JIRA_1521_Batch2-barcode06/intermediate/call_output/methylation_calls.cpg.tsv.gz
    [DEBUG]: Closing file:/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_1_date_8_3-20200908_P110_JIRA_1521_Batch1-barcode02/intermediate/call_output/methylation_calls.cpg.tsv.gz
    [DEBUG]: Closing file:/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_6_date_8_3-20200916_P111_JIRA_1521_Batch2-barcode08/intermediate/call_output/methylation_calls.cpg.tsv.gz
    [DEBUG]: Closing file:/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_1_date_8_3-20200918_P115_JIRA_1521_Batch3_reload-barcode09/intermediate/call_output/methylation_calls.cpg.tsv.gz
    [DEBUG]: Closing file:/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_3_date_8_3-20200918_P115_JIRA_1521_Batch3_reload-barcode10/intermediate/call_output/methylation_calls.cpg.tsv.gz
    [DEBUG]: Closing file:/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_1_date_8_3-20201207_P125_JIRA_1521_finalbatch_shear-barcode20/intermediate/call_output/methylation_calls.cpg.tsv.gz
    [DEBUG]: Closing file:/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_1_date_8_3-20201023_P123_JIRA_1521_Batch10_sheared-barcode20/intermediate/call_output/methylation_calls.cpg.tsv.gz
    [DEBUG]: Closing file:/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_5_date_8_3-20200928_P117_JIRA_1521_Batch5-barcode18/intermediate/call_output/methylation_calls.cpg.tsv.gz
    [DEBUG]: Closing file:/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_2_date_8_3-20201008_P118_JIRA_1521_Batch6_sheared-barcode02/intermediate/call_output/methylation_calls.cpg.tsv.gz
    [DEBUG]: Closing file:/local/ifs2_projdata/8570/projects/EPIENV3/jespinoz/Analysis/SOLEXASEQ-1521/f5c_output/rat_5_date_8_3-20201207_P125_JIRA_1521_finalbatch_shear-barcode18/intermediate/call_output/methylation_calls.cpg.tsv.gz
Traceback (most recent call last):
  File "/usr/local/devel/ANNOTATION/jespinoz/anaconda3/envs/ont_env/bin/pycoMeth", line 8, in <module>
    sys.exit(main())
  File "/usr/local/devel/ANNOTATION/jespinoz/anaconda3/envs/ont_env/lib/python3.7/site-packages/pycoMeth/__main__.py", line 116, in main
    args.func(**vars(args))
  File "/usr/local/devel/ANNOTATION/jespinoz/anaconda3/envs/ont_env/lib/python3.7/site-packages/pycoMeth/CpG_Aggregate.py", line 69, in CpG_Aggregate
    raise pycoMethError("Invalid input file type passed (nanopolish_fn). Expecting Nanopolish call_methylation output TSV file")
pycoMeth.common.pycoMethError: Invalid input file type passed (nanopolish_fn). Expecting Nanopolish call_methylation output TSV file
a-slide commented 3 years ago

Hi @jolespin pycoMeth is expecting the following field from nanopolish call-methylation output: "chromosome", "strand", "start", "end" , "read_name", "log_lik_ratio", "log_lik_methylated", "log_lik_unmethylated" Looking at the log, it doesn't seem like you have a "strand" field. What version of Nanopolish do you use ?

a-slide commented 3 years ago

Gzip format is not an issue. You can see that the files are actually all opened in gzip mode

jolespin commented 3 years ago

I ended up using the optimized reimplementation:

https://github.com/hasindu2008/f5c

It’s supposed to be a drop in replacement but I’ll check to see which fields are not consistent.

@hasindu2008 do you have any insight?

I could parse the bam files to add this column manually maybe.

On Jan 23, 2021, at 6:52 AM, Adrien Leger notifications@github.com wrote:

 Gzip format is not an issue. You can see that the files are actually all opened in gzip mode

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

hasindu2008 commented 3 years ago

@jolespin when you call methylation add --meth-out-version 2 to f5c to get the same fields as the latest nanopolish versions.

--meth-out-version INT Format version of the output Methylation tsv file If set to 1, the columns printed adhere to the output format of Nanopolish early versions. If set to 2, adhere to the latest nanopolish output format that additionally includes the strand column and the header num_cpgs renamed to num_motifs) [default value: 1]

jolespin commented 3 years ago

@jolespin

when you call methylation add --meth-out-version 2 to f5c to get the same fields as the latest nanopolish versions.

--meth-out-version INT

Format version of the output Methylation tsv file

If set to 1, the columns printed adhere to the output format of Nanopolish early versions. If set to 2, adhere to the latest nanopolish output format that additionally includes the strand column and the header num_cpgs renamed to num_motifs) [default value: 1]

Hmm... I just ran f5c on hundreds of samples. Is there any way I can convert?

Alternatively, is there an earlier version of pycometh that handles earlier versions of Nanopolish?

hasindu2008 commented 3 years ago

@jolespin An year ago, I ran pycoMeth (pycoMeth v0.3.4) on such samples by manually renaming the num_cpgs to num_motifs in the f5c tsv outputs. Not sure if the strand column is used in latest pycometh for any calculations, but at that time it was not required I guess.

After renaming they looked like:

chromosome      start   end     read_name       log_lik_ratio   log_lik_methylated      log_lik_unmethylated    num_calling_strands     num_motifs      sequence
chr2    124998  124998  72d14342-4eb2-4d6e-8a8e-bd65ea179dc2    4.26    -111.34 -115.60 1       1       ACCTGCGAACA

Command was pycoMeth CpG_Aggregate -i $input -f hg38noAlt.fa -t $each.tsv -b $each.bed --progress

jolespin commented 3 years ago

@hasindu2008 Thanks, I got it to run w/ pycoMeth v03.4 after relabeling num_cpgs to num_motifs. I've changed my f5c wrapper to use --meth-out-version 2 by default to avoid this in the future.

Do you know of a quicker way to relabel the header? I did a very naive way for my first batch of files by loading everything into pandas (Python), relabeling, and then writing to disk.