Nicole's moPepGen adventure #746

alkaZeltser · 2023-05-19T18:56:07Z

alkaZeltser
May 19, 2023
Collaborator

I'll keep track of my progress with test-running moPepGen here.

Starting inputs: Roni's radioresistance cell line dataset /hot/project/disease/ProstateTumor/PRAD-000096-RadioResDU145Molecular/

alkaZeltser · 2023-05-19T19:07:23Z

alkaZeltser
May 19, 2023
Collaborator Author

Installation

I started by opening up my favorite miniconda env and running
pip install git+ssh://[email protected]:uclahs-cds/private-moPepGen.git

This did not work with the following error:

Collecting git+ssh://****@github.com:uclahs-cds/private-moPepGen.git
  Cloning ssh://****@github.com:uclahs-cds/private-moPepGen.git to /tmp/pip-req-build-_mdrhz40
  Running command git clone -q 'ssh://****@github.com:uclahs-cds/private-moPepGen.git' /tmp/pip-req-build-_mdrhz40
  ssh: Could not resolve hostname github.com:uclahs-cds: Name or service not known
  fatal: Could not read from remote repository.

  Please make sure you have the correct access rights
  and the repository exists.

WARNING: Discarding git+ssh://****@github.com:uclahs-cds/private-moPepGen.git. Command errored out with exit status 128: git clone -q 'ssh://****@github.com:uclahs-cds/private-moPepGen.git' /tmp/pip-req-build-_mdrhz40 Check the logs for full command output.
ERROR: Command errored out with exit status 128: git clone -q 'ssh://****@github.com:uclahs-cds/private-moPepGen.git' /tmp/pip-req-build-_mdrhz40 Check the logs for full command output.

I then tried to just clone the repo, which worked with no issues:
git clone [email protected]:uclahs-cds/private-moPepGen.git

So I have all the code now but I don't know how to actually run it from command line - it's not officially installed and I don't know what binary I'm supposed to be executing.

5 replies

zhuchcn May 20, 2023
Maintainer

There is an issue with the url. Try this command:

pip install git+ssh://[email protected]/uclahs-cds/private-moPepGen.git

The url style requried by pip is different from the ssh remote url on github and I alway got confused. Btw, will be good to install it in a new conda environment so the dependency versions won't have conflict. Alternatively, we have the docker image available, so just do this:

docker pull ghcr.io/uclahs-cds/mopepgen:dev

alkaZeltser May 20, 2023
Collaborator Author

I started a fresh conda env and installed using
pip install git+ssh://[email protected]/uclahs-cds/private-moPepGen.git

which successfully built moPepGen and installed dependencies. It took longer than I expected, I think close to 10 minutes - not sure I lost track. But successful nonetheless.

(moPep_env) [nzeltser@ip-0A125209 moPepGen]$ moPepGen
usage: moPopGen [-h] [-V] <command> [options]

optional arguments:
  -h, --help     show this help message and exit
  -V, --version  Show version number and exit.

subcommands:
  
  -- Indexing
     generateIndex       Generate genome and proteome index files for moPepGen.
     indexGVF            Generate index for GVF files.
  
  -- Parsing
     parseVEP            Parse VEP output.
     parseREDItools      Parse REDItools annotated output.
     parseSTARFusion     Parse STAR-Fusion output.
     parseFusionCatcher  Parse FusionCatcher output.
     parseArriba         Parse Arriba output.
     parseRMATS          Parse rMATS output.
     parseCIRCexplorer   Parse CIRCexplorer known circRNA output.
  
  -- Calling
     callVariant         Call non-canonical peptides from genomic variants.
     callNoncoding       Call non-canonical peptides from noncoding transcripts.
     callAltTranslation  Call non-canonital peptides with alternative translation
                         from coding transcripts.
  
  -- Processing
     summarizeFasta      Summarize variant peptides.
     filterFasta         Filter noncanonical peptides.
     splitFasta          Split noncanonical peptides into separate databases.
     mergeFasta          Merge multiple variant peptide FASTA databases.
     encodeFasta         Encode variant peptide FASTA file header.
     decoyFasta          Generate decoy database FASTA file.

On to the next step!

zhuchcn May 21, 2023
Maintainer

Do you remember which step took 10 min? I'm guessing it's either setting up conda, or downloading it because of network speed. Installing it with pip should be pretty quick

alkaZeltser May 21, 2023
Collaborator Author

It was the entire pip install command. I'm not counting the conda setup. The command didn't return for at least 10 minutes.

zhuchcn May 30, 2023
Maintainer

It used to install pretty fast before we added the dependency of matplotlib, which also depends on a bunch of heavy packages such as numpy. So I guess the majority of time was spent on installing matplotlib & numpy.

alkaZeltser · 2023-05-21T00:07:24Z

alkaZeltser
May 21, 2023
Collaborator Author

Reference Index

The instructions are very clear that the first step is to generate an index of reference files for moPepGen. It's not clear whether this reference set needs to be built from the same set of reference files that was used to call your genomic inputs. If that is the case, then it might be important to note that DNA alignment and RNA processing was not necessarily performed on the same reference set. In my case DNA alignment was on GRCh38-BI-20160721 and RNA was on GRCh38-EBI-GENCODE36. I don't actually know the difference between the BWA indexed reference and the standard GENCODE36.

For now, pushing on with the GENCODE reference set:

moPepGen generateIndex \
-g /hot/ref/reference/GRCh38-EBI-GENCODE36/GRCh38.p13.genome.fa \
-a /hot/ref/reference/GRCh38-EBI-GENCODE36/gencode.v36.chr_patch_hapl_scaff.annotation.gtf \
-p /hot/ref/reference/GRCh38-EBI-GENCODE36/gencode.v36.pc_translations.fa \
-o ./index

Initially tried this in an F2 node, got a memory error, fair enough. Might be nice to give an idea of the mem requirements.

[ 2023-05-20 16:39:02 ] moPepGen generateIndex started
Traceback (most recent call last):
  File "/hot/users/nzeltser/miniconda3/envs/moPep_env/bin/moPepGen", line 8, in <module>
    sys.exit(main())
  File "/hot/users/nzeltser/miniconda3/envs/moPep_env/lib/python3.9/site-packages/moPepGen/cli/__main__.py", line 92, in main
    args.func(args)
  File "/hot/users/nzeltser/miniconda3/envs/moPep_env/lib/python3.9/site-packages/moPepGen/cli/generate_index.py", line 71, in generate_index
    genome.dump_fasta(path_genome)
  File "/hot/users/nzeltser/miniconda3/envs/moPep_env/lib/python3.9/site-packages/moPepGen/dna/DNASeqDict.py", line 37, in dump_fasta
    for record in SeqIO.parse(path, 'fasta'):
  File "/hot/users/nzeltser/miniconda3/envs/moPep_env/lib/python3.9/site-packages/Bio/SeqIO/Interfaces.py", line 72, in __next__
    return next(self.records)
  File "/hot/users/nzeltser/miniconda3/envs/moPep_env/lib/python3.9/site-packages/Bio/SeqIO/FastaIO.py", line 238, in iterate
    for title, sequence in SimpleFastaParser(handle):
  File "/hot/users/nzeltser/miniconda3/envs/moPep_env/lib/python3.9/site-packages/Bio/SeqIO/FastaIO.py", line 65, in SimpleFastaParser
    yield title, "".join(lines).replace(" ", "").replace("\r", "")
MemoryError

Tried again in an f16, got a little farther, then it crashed without much explanation. Note the warning from Biopython:

[ 2023-05-20 16:44:07 ] moPepGen generateIndex started
[ 2023-05-20 16:44:44 ] Genome FASTA loaded
[ 2023-05-20 16:45:09 ] Genome FASTA saved to disk.
/hot/users/nzeltser/miniconda3/envs/moPep_env/lib/python3.9/site-packages/Bio/SeqFeature.py:257: BiopythonDeprecationWarning: Using the strand argument is deprecated, and will be removed in a future release. Please set it via the location object instead.
  warnings.warn(
Killed

Next I ramped it up to an F72:

sbatch --exclusive --partition=F72 -J moIndex -e moIndex.err -o moIndex.log --wrap="moPepGen generateIndex -g /hot/ref/reference/GRCh38-EBI-GENCODE36/GRCh38.p13.genome.fa -a /hot/ref/reference/GRCh38-EBI-GENCODE36/gencode.v36.chr_patch_hapl_scaff.annotation.gtf -p /hot/ref/reference/GRCh38-EBI-GENCODE36/gencode.v36.pc_translations.fa -o ./index"

Success!

Log:

[ 2023-05-20 16:55:57 ] moPepGen generateIndex started
[ 2023-05-20 16:56:33 ] Genome FASTA loaded
[ 2023-05-20 16:56:48 ] Genome FASTA saved to disk.
[ 2023-05-20 16:58:51 ] Genome annotation GTF loaded.
[ 2023-05-20 17:00:13 ] Genome annotation GTF saved to disk.
[ 2023-05-20 17:00:13 ] Proteome FASTA loaded.
[ 2023-05-20 17:00:14 ] Proteome FASTA saved to disk.
[ 2023-05-20 17:02:27 ] canonical peptide pool generated.
[ 2023-05-20 17:02:28 ] canonical peptide pool saved to disk.
[ 2023-05-20 17:02:32 ] rss=230.7 MiB, vms=2.847 GiB, wallclock=0:05:48.790000, system=0:00:19.640000, cpu_usage=91.5%

Err:

/hot/users/nzeltser/miniconda3/envs/moPep_env/lib/python3.9/site-packages/Bio/SeqFeature.py:257: BiopythonDeprecationWarning: Using the strand argument is deprecated, and will be removed in a future release. Please set it via the location object instead.
  warnings.warn(

Outputs:

(moPep_env) [nzeltser@ip-0A12525A moPepGen]$ ls -lh index/
total 3.2G
-rw-r--r--. 1 nzeltser domain users  18M May 20 17:00 annotation.dat
-rw-r--r--. 1 nzeltser domain users  51M May 20 17:02 canonical_peptides.pkl
-rw-r--r--. 1 nzeltser domain users 2.0M May 20 17:02 coding_transcripts.pkl
-rw-r--r--. 1 nzeltser domain users 3.1G May 20 16:56 genome.pkl
-rw-r--r--. 1 nzeltser domain users  64M May 20 17:00 proteome.pkl

tl;dr

would like more clarity on relationship between inputs and references
Biopython Depcreation Warning
might be nice to provide an idea of memory requirements

1 reply

zhuchcn May 21, 2023
Maintainer

In my case DNA alignment was on GRCh38-BI-20160721 and RNA was on GRCh38-EBI-GENCODE36. I don't actually know the difference between the BWA indexed reference and the standard GENCODE36.

For DNA alignment, it's fine as long as the reference genome's major version matches (GRCh38). But when you run VEP, make sure you use the GTF file from GENCODE36.

Biopython Depcreation Warning

We should also evaluate the dependencies. I sometimes do see different python/biopython version may cause different issues. But using our docker image is guaranteed to be fine!

might be nice to provide an idea of memory requirements

That's a good point. We should evaluate all commands and document the approximate memorry usage of each.

alkaZeltser · 2023-05-21T00:28:03Z

alkaZeltser
May 21, 2023
Collaborator Author

Question re input file format

For all moPepGen functions, can -i receive a file that is a list of paths, or does each input file have to be one path in one line?

1 reply

zhuchcn May 21, 2023
Maintainer

It depends on the command. callVariant and parseVEP are the only two commands that their -i can accept multiple files. It's easier to look at the usage. If its like this, it accepts multiple (callVariant)

-i ['<files>'] [['<files>'] ...], --input-path ['<files>'] [['<files>'] ...]

And this is when the command only accept one file (parseSTARFusion).

-i <file>, --input-path <file>

alkaZeltser · 2023-05-29T07:14:59Z

alkaZeltser
May 29, 2023
Collaborator Author

Parsing Part 1.1 - Fusion from STAR-Fusion

I'm starting the GVF file creation. I wanted to start with SNVs but turns out our annotation pipeline doesn't actually have VEP, despite documentation saying otherwise LOL. So starting with something easier for now: Star fusion. No problems.

Code for Star Fusion parsing for all available samples:

INPUT_LIST="/hot/user/nzeltser/moPepGen/Fusion/paths-to-star-fusion-results.txt"
OUTPUT_DIR="/hot/user/nzeltser/moPepGen/Fusion/output"
INDEX_DIR="/hot/user/nzeltser/moPepGen/index"

# loop through paths in INPUT_LIST text file and run a STAR-Fusion parse on each
while read -r line; do
    echo "Parsing STAR-Fusion output for $line"
    # parse file basename without extension
    filename=$(basename -- "$line")
    moPepGen parseSTARFusion \
        -i $line \
        --index-dir $INDEX_DIR \
        --source Fusion \
        -o ${OUTPUT_DIR}/${filename}_fusion.gvf
done < "$INPUT_LIST"

Sample from log:

Parsing STAR-Fusion output for /hot/project/disease/ProstateTumor/PRAD-000096-RadioResDU145Molecular/fusion_trancript/star_results/LHPRCRRP000001-T001-C01-F_star-fusion.fusion_predictions.tsv
[ 2023-05-28 22:11:31 ] moPepGen parseSTARFusion started
[ 2023-05-28 22:13:26 ] Reference indices loaded.
[ 2023-05-28 22:13:26 ] STAR-Fusion output /hot/project/disease/ProstateTumor/PRAD-000096-RadioResDU145Molecular/fusion_trancript/star_results/LHPRCRRP000001-T001-C01-F_star-fusion.fusion_predictions.tsv loaded.
[ 2023-05-28 22:13:26 ] Variants sorted.
[ 2023-05-28 22:13:26 ] Variant info written to disk.
[ 2023-05-28 22:13:28 ] rss=187 MiB, vms=2.801 GiB, wallclock=0:01:38.320000, system=0:00:08.950000, cpu_usage=85.4%

Sample from err:

/hot/users/nzeltser/miniconda3/envs/moPep_env/lib/python3.9/site-packages/Bio/SeqFeature.py:257: BiopythonDeprecationWarning: Using the strand argument is deprecated, and will be removed in a future release. Please set it via the location object instead.
  warnings.warn(

^Biopython Deprecation warning again.

3 replies

zhuchcn May 29, 2023
Maintainer

Yup the SeqFeature strand issue is known. SeqFeature has a location:FeatureLocation attribute, which has strand already, so I think biopython wants to avoid strand having different values from SeqFeature.strand and SeqFeature.lcation.strand.

zhuchcn May 29, 2023
Maintainer

We have an issue opened at #616

zhuchcn May 30, 2023
Maintainer

This should be fixed in #752

alkaZeltser · 2023-05-30T23:33:50Z

alkaZeltser
May 30, 2023
Collaborator Author

Parsing Part 1.2 - Fusion from Fusion Catcher

This ran very smoothly, no issues.

Parsing code:

INPUT_LIST="/hot/user/nzeltser/moPepGen/Fusion-catcher/paths-to-fusion-catcher-results.txt"
OUTPUT_DIR="/hot/user/nzeltser/moPepGen/Fusion-catcher/output"
INDEX_DIR="/hot/user/nzeltser/moPepGen/index"

# loop through paths in INPUT_LIST text file and run a STAR-Fusion parse on each
while read -r line; do
    echo "Parsing fusion-catcher output for $line"
    # parse file basename without extension
    filename=$(basename -- "$line")
    moPepGen parseFusionCatcher \
        -i $line \
        --index-dir $INDEX_DIR \
        --source Fusion \
        -o ${OUTPUT_DIR}/${filename}_fusion.gvf
done < "$INPUT_LIST"

Log:

Parsing fusion-catcher output for /hot/project/disease/ProstateTumor/PRAD-000096-RadioResDU145Molecular/fusion_trancript/fusioncatcher_results/LHPRCRRP000001-T001-C01-F_final-list_candidate-fusion-genes.txt
[ 2023-05-30 15:48:13 ] moPepGen parseFusionCatcher started
[ 2023-05-30 15:51:44 ] Reference indices loaded.
[ 2023-05-30 15:51:44 ] FusionCatcher output /hot/project/disease/ProstateTumor/PRAD-000096-RadioResDU145Molecular/fusion_trancript/fusioncatcher_results/LHPRCRRP000001-T001-C01-F_final-list_candidate-fusion-genes.txt loaded.
[ 2023-05-30 15:51:44 ] Variants sorted.
[ 2023-05-30 15:51:44 ] Variants written to disk.
[ 2023-05-30 15:51:46 ] rss=188.1 MiB, vms=2.804 GiB, wallclock=0:01:47.780000, system=0:00:11.460000, cpu_usage=52.4%

Error (being addressed in #752) :

/hot/users/nzeltser/miniconda3/envs/moPep_env/lib/python3.9/site-packages/Bio/SeqFeature.py:257: BiopythonDeprecationWarning: Using the strand argument is deprecated, and will be removed in a future release. Please set it via the location object instead.
  warnings.warn(

0 replies

alkaZeltser · 2023-05-31T01:03:57Z

alkaZeltser
May 31, 2023
Collaborator Author

Parsing Part 3 - Alternative Splicing

Testing functionality of moPepGen parseRMATS. No problems.

Code:

while read -r line; do
    echo "Parsing rMATS output for $line"
    # Parse file basename without extension and properly interpret input path with wildcards
    filename=$(find "$line" -path "$line/2021*/rmats-v4.1.1/*SE.MATS.JC.txt" -exec basename {} \; | cut -d_ -f1)
    moPepGen parseRMATS \
        --se $line/2021*/rmats-v4.1.1/*SE.MATS.JC.txt \
        --a3ss $line/2021*/rmats-v4.1.1/*A3SS.MATS.JC.txt \
        --a5ss $line/2021*/rmats-v4.1.1/*A5SS.MATS.JC.txt \
        --mxe $line/2021*/rmats-v4.1.1/*MXE.MATS.JC.txt \
        --ri $line/2021*/rmats-v4.1.1/*RI.MATS.JC.txt \
        --index-dir $INDEX_DIR \
        -o ${OUTPUT_DIR}/${filename}_alt_splice.gvf \
        --source AltSplice
done < "$INPUT_LIST"

Logs:

Parsing rMATS output for /hot/project/disease/ProstateTumor/PRAD-000096-RadioResDU145Molecular/SpliceIsoform/LHPRCRRP000003-T003-C03-F
[ 2023-05-30 17:55:40 ] moPepGen parseRMATS started
[ 2023-05-30 17:57:19 ] Reference indices loaded.
[ 2023-05-30 17:57:19 ] Start parsing SE file /hot/project/disease/ProstateTumor/PRAD-000096-RadioResDU145Molecular/SpliceIsoform/LHPRCRRP000003-T003-C03-F/20210813-183746/rmats-v4.1.1/LHPRCRRP000003-T003-C03-F_SE.MATS.JC.txt
[ 2023-05-30 17:57:29 ] Start parsing A5SS file /hot/project/disease/ProstateTumor/PRAD-000096-RadioResDU145Molecular/SpliceIsoform/LHPRCRRP000003-T003-C03-F/20210813-183746/rmats-v4.1.1/LHPRCRRP000003-T003-C03-F_A5SS.MATS.JC.txt
[ 2023-05-30 17:57:29 ] Start parsing A3SS file /hot/project/disease/ProstateTumor/PRAD-000096-RadioResDU145Molecular/SpliceIsoform/LHPRCRRP000003-T003-C03-F/20210813-183746/rmats-v4.1.1/LHPRCRRP000003-T003-C03-F_A3SS.MATS.JC.txt
[ 2023-05-30 17:57:31 ] Start parsing MXE file /hot/project/disease/ProstateTumor/PRAD-000096-RadioResDU145Molecular/SpliceIsoform/LHPRCRRP000003-T003-C03-F/20210813-183746/rmats-v4.1.1/LHPRCRRP000003-T003-C03-F_MXE.MATS.JC.txt
[ 2023-05-30 17:57:32 ] Start parsing RI file /hot/project/disease/ProstateTumor/PRAD-000096-RadioResDU145Molecular/SpliceIsoform/LHPRCRRP000003-T003-C03-F/20210813-183746/rmats-v4.1.1/LHPRCRRP000003-T003-C03-F_RI.MATS.JC.txt
[ 2023-05-30 17:57:34 ] Variants sorted.
[ 2023-05-30 17:57:34 ] Variants written to disk.
[ 2023-05-30 17:57:37 ] rss=190 MiB, vms=2.814 GiB, wallclock=0:01:54.340000, system=0:00:09.420000, cpu_usage=99.9%

Error, same as before:

/hot/users/nzeltser/miniconda3/envs/moPep_env/lib/python3.9/site-packages/Bio/SeqFeature.py:257: BiopythonDeprecationWarning: Using the strand argument is deprecated, and will be removed in a future release. Please set it via the location object instead.
  warnings.warn(

1 reply

zhuchcn Jun 23, 2023
Maintainer

Error, same as before:

What docker version are you using? I thought I fixed it but maybe this is a different issue.

alkaZeltser · 2023-06-21T14:44:55Z

alkaZeltser
Jun 21, 2023
Collaborator Author

Parsing Part 4: SNVs + VEP

After many tribulations with VEP, I think I have finally achieved a VEP parse.

VEP Annotation

Lessons learned

GTF file provided to VEP must be sorted, bgzipped, and indexed
If specifying chromosomes via --chr, must make sure to match chr name format to VCF e.g. "chr1" instead of "1"
The Boutros lab VEP docker was producing an error that I couldn't resolve, I switched to the ensembl docker which didn't produce an error

Code:

while read -r line; do
    echo "Running VEP on $line"
    # split file path from $line on slashes and extract seventh element
    filename=$(echo "$line" | cut -d'/' -f8)
    # run VEP
    docker run -u $(id -u):$(id -g) --rm -v $line:/in.vcf -v $OUTPUT_DIR:/out -v /hot/ref:/hot/ref \
    -v /hot/user/nzeltser/moPepGen/annotate-vep/inputs/ref:/hot/user/nzeltser/moPepGen/annotate-vep/inputs/ref \
    -v ${CACHE_DIR}:/.vep \
    ensemblorg/ensembl-vep \
        vep \
        --offline \
        --cache \
        --cache_version 104 \
        --fork ${N_THREADS} \
        --buffer_size 5000 \
        -i /in.vcf \
        -o /out/${filename}_vep.tsv \
        --assembly GRCh38 \
        --fasta ${FASTA} \
        --no_stats \
        --no_intergenic \
        --custom ${ANNOTATION_GTF},${REFERENCE_VERSION},gtf \
        --check_ref

    # # run vep filter
    echo "Running VEP filter on $line"
    docker run -u $(id -u):$(id -g) --rm -v $line:/in.vcf -v $OUTPUT_DIR:/out -v /hot/ref:/hot/ref \
    -v /hot/user/nzeltser/moPepGen/annotate-vep/inputs/ref:/hot/user/nzeltser/moPepGen/annotate-vep/inputs/ref \
    ghcr.io/uclahs-cds/vep:101.0 \
        filter_vep \
        --force_overwrite \
        -i /out/${filename}_vep.tsv \
        -o /out/${filename}_vep_filtered.tsv \
        --filter "Source = ${REFERENCE_VERSION}"

done < "$INPUT_LIST"

VEP Parsing:

No issues other than ye old Biopython warning /hot/users/nzeltser/miniconda3/envs/moPep_env/lib/python3.9/site-packages/Bio/SeqFeature.py:257: BiopythonDeprecationWarning: Using the strand argument is deprecated, and will be removed in a future release. Please set it via the location object instead.

Code:

while read -r line; do
    echo "Parsing VEP output for $line"
    # parse file basename without extension
    filename=$(basename -- "$line" | cut -d'_' -f1)
    moPepGen parseVEP \
        -i $line \
        --index-dir $INDEX_DIR \
        --source SNV \
        -o ${OUTPUT_DIR}/${filename}_fusion.gvf
done < "$INPUT_LIST"

Logs:

[ 2023-06-21 04:51:51 ] moPepGen parseVEP started
[ 2023-06-21 04:53:31 ] Reference indices loaded.
[ 2023-06-21 05:49:29 ] VEP file /hot/user/nzeltser/moPepGen/annotate-vep/output/VEP-ensembl/LHPRCRRP000003-T002-C02-F-LHPRCRRP000002-T006-C06-F_vep_filtered.tsv loaded.
[ 2023-06-21 05:49:43 ] VEP sorting done.
[ 2023-06-21 05:50:27 ] Variant info written to disk.
[ 2023-06-21 05:50:35 ] rss=304.3 MiB, vms=1.675 GiB, wallclock=0:57:20.920000, system=0:00:14.330000, cpu_usage=98.0%

2 replies

zhuchcn Jun 23, 2023
Maintainer

GTF file provided to VEP must be sorted, bgzipped, and indexed

If specifying chromosomes via --chr, must make sure to match chr name format to VCF e.g. "chr1" instead of "1"

Will document this.

The Boutros lab VEP docker was producing an error that I couldn't resolve, I switched to the ensembl docker which didn't produce an error

Do you still have the error message? I had the opposite thing. I had issue with ensembl's docker so I decided to create our own. Maybe I can fix that issue.

alkaZeltser Jul 3, 2023
Collaborator Author

This is the full error:

-------------------- EXCEPTION -------------------- 
MSG: ERROR: Forked process(es) died: read-through of cross-process communication detected 
STACK Bio::EnsEMBL::VEP::Runner::_forked_buffer_to_output /usr/local/share/ensembl-vep-101.0-1/modules/Bio/EnsEMBL/VEP/Runner.pm:565
STACK Bio::EnsEMBL::VEP::Runner::next_output_line /usr/local/share/ensembl-vep-101.0-1/modules/Bio/EnsEMBL/VEP/Runner.pm:370
STACK Bio::EnsEMBL::VEP::Runner::run /usr/local/share/ensembl-vep-101.0-1/modules/Bio/EnsEMBL/VEP/Runner.pm:211
STACK toplevel /usr/local/bin/vep:228
Date (localtime) = Thu Jun 15 23:27:41 2023 Ensembl API version = 101
---------------------------------------------------
Possible precedence issue with control flow operator at /usr/local/lib/site_perl/5.26.2/Bio/DB/IndexedBase.pm line 805.

alkaZeltser · 2023-07-03T22:36:51Z

alkaZeltser
Jul 3, 2023
Collaborator Author

Question about GTF reference files:

So I finished the VEP annotation and am moving onto RNA Editing Sites annotation, and both of these require a reference annotation GTF. I just realized that the GTF I used for VEP annotation had to be sorted and bgzipped. However the GTF used in MoPep indexing was not a sorted or bgzipped version:

moPepGen generateIndex \
-g /hot/ref/reference/GRCh38-EBI-GENCODE36/GRCh38.p13.genome.fa \
-a /hot/ref/reference/GRCh38-EBI-GENCODE36/gencode.v36.chr_patch_hapl_scaff.annotation.gtf \ # THIS GTF IS NOT SORTED
-p /hot/ref/reference/GRCh38-EBI-GENCODE36/gencode.v36.pc_translations.fa \
-o ./index

vep \
        --offline \
        --cache \
        --cache_version 104 \
        --fork ${N_THREADS} \
        --buffer_size 5000 \
        -i /in.vcf \
        -o /out/${filename}_vep.tsv \
        --assembly GRCh38 \
        --fasta ${FASTA} \
        --no_stats \
        --no_intergenic \
        --custom ${ANNOTATION_GTF},${REFERENCE_VERSION},gtf \ # THIS GTF IS SORTED
        --check_ref

Is that going to be an issue? Does the indexed GTF need to be completely identical, even down to the sorting of the contents?
If so, then I might have to re-do the indexing step and all of my other parsing steps (they all take the --index-dir argument, and that dir would be containing a non-matching GTF).

5 replies

zhuchcn Jul 3, 2023
Maintainer

The GTF file being passed to generateIndex should not be sorted or bgzipped. Just sorting it may not be an issue but bgzip is definitely going to be.

Btw, it may not be a bad idea to redo the indexing step with the new 1.1.0, in which we significantly reduced memory usage, even for parsers.

alkaZeltser Jul 3, 2023
Collaborator Author

Err okay I'm very confused, you threw in a double negative wrapping an "or" statement LOL.

So, the GTF file passed to generateIndex needs to be

sorted
not bgzipped
???
BTW the link in this line in the Vignette is broken:
https://probable-adventure-7wqgn3w.pages.github.io/vignette/#:~:text=See%20here%20for%20more%20details.

alkaZeltser Jul 3, 2023
Collaborator Author

Also, I'm not using the MoPepGen docker, I installed it into a conda env per Vignette instructions. Would you like to add instructions for updating that installation to the Vignette?

zhuchcn Jul 4, 2023
Maintainer

Just use the GTF directly from gencode/ensembl!

Which line is broken? Updating is the same as installing. So just do this:

pip install git+ssh://[email protected]/uclahs-cds/[email protected]

zhuchcn Jul 4, 2023
Maintainer

"should not be sorted or bgzipped" is what I meant to say. Updated already

alkaZeltser · 2023-07-04T01:59:27Z

alkaZeltser
Jul 4, 2023
Collaborator Author

In the Reference Data section, the last line in the paragraph has a hyperlink: "See here for more details. " The link in "here" is broken, leads to 404 not found.

2 replies

zhuchcn Jul 5, 2023
Maintainer

This should be fixed in #763. Can you take a look at the PR?

alkaZeltser Jul 5, 2023
Collaborator Author

Yep, looks good

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Nicole's moPepGen adventure #746

{{title}}

Replies: 9 comments 20 replies

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

Select a reply

Nicole's moPepGen adventure #746

alkaZeltser May 19, 2023 Collaborator

Replies: 9 comments · 20 replies

alkaZeltser May 19, 2023 Collaborator Author

Installation

zhuchcn May 20, 2023 Maintainer

alkaZeltser May 20, 2023 Collaborator Author

zhuchcn May 21, 2023 Maintainer

alkaZeltser May 21, 2023 Collaborator Author

zhuchcn May 30, 2023 Maintainer

alkaZeltser May 21, 2023 Collaborator Author

Reference Index

tl;dr

zhuchcn May 21, 2023 Maintainer

alkaZeltser May 21, 2023 Collaborator Author

Question re input file format

zhuchcn May 21, 2023 Maintainer

alkaZeltser May 29, 2023 Collaborator Author

Parsing Part 1.1 - Fusion from STAR-Fusion

zhuchcn May 29, 2023 Maintainer

zhuchcn May 29, 2023 Maintainer

zhuchcn May 30, 2023 Maintainer

alkaZeltser May 30, 2023 Collaborator Author

Parsing Part 1.2 - Fusion from Fusion Catcher

alkaZeltser May 31, 2023 Collaborator Author

Parsing Part 3 - Alternative Splicing

zhuchcn Jun 23, 2023 Maintainer

alkaZeltser Jun 21, 2023 Collaborator Author

Parsing Part 4: SNVs + VEP

VEP Annotation

VEP Parsing:

zhuchcn Jun 23, 2023 Maintainer

alkaZeltser Jul 3, 2023 Collaborator Author

alkaZeltser Jul 3, 2023 Collaborator Author

zhuchcn Jul 3, 2023 Maintainer

alkaZeltser Jul 3, 2023 Collaborator Author

alkaZeltser Jul 3, 2023 Collaborator Author

zhuchcn Jul 4, 2023 Maintainer

zhuchcn Jul 4, 2023 Maintainer

alkaZeltser Jul 4, 2023 Collaborator Author

zhuchcn Jul 5, 2023 Maintainer

alkaZeltser Jul 5, 2023 Collaborator Author

alkaZeltser
May 19, 2023
Collaborator

Replies: 9 comments 20 replies

alkaZeltser
May 19, 2023
Collaborator Author

zhuchcn May 20, 2023
Maintainer

alkaZeltser May 20, 2023
Collaborator Author

zhuchcn May 21, 2023
Maintainer

alkaZeltser May 21, 2023
Collaborator Author

zhuchcn May 30, 2023
Maintainer

alkaZeltser
May 21, 2023
Collaborator Author

zhuchcn May 21, 2023
Maintainer

alkaZeltser
May 21, 2023
Collaborator Author

zhuchcn May 21, 2023
Maintainer

alkaZeltser
May 29, 2023
Collaborator Author

zhuchcn May 29, 2023
Maintainer

zhuchcn May 29, 2023
Maintainer

zhuchcn May 30, 2023
Maintainer

alkaZeltser
May 30, 2023
Collaborator Author

alkaZeltser
May 31, 2023
Collaborator Author

zhuchcn Jun 23, 2023
Maintainer

alkaZeltser
Jun 21, 2023
Collaborator Author

zhuchcn Jun 23, 2023
Maintainer

alkaZeltser Jul 3, 2023
Collaborator Author

alkaZeltser
Jul 3, 2023
Collaborator Author

zhuchcn Jul 3, 2023
Maintainer

alkaZeltser Jul 3, 2023
Collaborator Author

alkaZeltser Jul 3, 2023
Collaborator Author

zhuchcn Jul 4, 2023
Maintainer

zhuchcn Jul 4, 2023
Maintainer

alkaZeltser
Jul 4, 2023
Collaborator Author

zhuchcn Jul 5, 2023
Maintainer

alkaZeltser Jul 5, 2023
Collaborator Author