chengl7-lab / scape

A package for estimating alternative polyadenylation events from scRNA-seq data.
MIT License
5 stars 1 forks source link

Getting empty files when running scape prepare_input #3

Closed msaland closed 4 months ago

msaland commented 4 months ago

Hello,

I'm running _scape prepare_input --utr_file GRCm39_112/GRCm39_112.csv --cb_file ${file}-barcodes.tsv --bam_file ${file}-sorted.bam --outputdir ${file}/ --chunksize 100 and I'm just getting empty .pkl files and nothing being generated. It's generating the barcode_index.csv fine though.

This is the output messages I'm getting:

2024-07-03 21:04:40,552 - INFO - All pickle files respective to file-sorted in pkl_input are deleted. 2024-07-03 21:04:40,758 - INFO - file/pkl_input/file-sorted.100.tmp.1.input.pkl is in processed 2024-07-03 21:31:24,335 - INFO - FINISHED 2024-07-03 21:31:24,335 - INFO - There are 1 pickle input files for this BAM file file-sorted.bam 2024-07-03 21:31:24,335 - INFO - Phrase tmp in name of pickle files is replaced by 1

Processing file-sorted.bam Create index file for barcode file/barcode_index.csv Done 0 files in 28.85623295853335 minutes.

I'm not sure what the issue is, but here is what the first few lines of my bam file looks like:

CTCAGAATAAACTATGCAAT,ACATGGGG,LH00403:55:223CL5LT1:1:1150:31965:9444 163 1 3149714 255 1S28M1S = 3149714 28 GTGACAAGAAGAAGAAGAAGAAGATTAAGG IIIIIIIIIIIIIIIIIIIIIIIIIIIIII NH:i:1 HI:i:1 AS:i:47 nM:i:4 CB:Z:CTCAGAATAAACTATGCAAT UR:Z:ACATGGGG UY:Z:JJJJJJJJJJ UB:Z:ACATGGGG CTCAGAATAAACTATGCAAT,ACATGGGG,LH00403:55:223CL5LT1:1:1150:31965:9444 83 1 3149714 255 1S28M1S = 3149714 -28 GTGACAAGAAGAAGAAGAAGAAGATTAAGG IIIIIIIIIIIIIIIIIIIIII9IIIIIII NH:i:1 HI:i:1 AS:i:47 nM:i:4 CB:Z:CTCAGAATAAACTATGCAAT UR:Z:ACATGGGG UY:Z:JJJJJJJJJJ UB:Z:ACATGGGG CAGCCATGATTTACCGAGGC,GGTGTTTG,LH00403:55:223CL5LT1:1:1216:31253:3532 163 1 3155985 255 135M = 3289229 133262 GGTAATAGGAGTAGATGAGAATGAAGATTTTCAACTTAAAGGGCCAGCAAATATCTTCAACAAAATTATAGAAGAAAACTTCCCAAACCTAAAGAAAGAGATGCCCATGAAAATACAAGAAGCCTACAGAACTCC IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII-IIIIIIIIIII NH:i:1 HI:i:1 AS:i:139 nM:i:5 CB:Z:CAGCCATGATTTACCGAGGC UR:Z:GGTGTTTG UY:Z:JJJJJJJJJJ UB:Z:GGTGTTTG TTCTGGCGCAACTATGCAAT,TGACTTAT,LH00403:55:223CL5LT1:1:2119:38484:25258 163 1 3172276 255 29M2S = 3172276 29 GATACAGAGAAGATTAGCATGGGCCCTGCGC -IIIIIIIIIIIIIIIIIIIIIIIIIIIIII NH:i:1 HI:i:1 AS:i:57 nM:i:0 CB:Z:TTCTGGCGCAACTATGCAAT UR:Z:TGACTTAT UY:Z:JJJJJJJJJJ UB:Z:TGACTTAT TTCTGGCGCAACTATGCAAT,TGACTTAT,LH00403:55:223CL5LT1:1:2119:38484:25258 83 1 3172276 255 29M2S = 3172276 -29 GATACAGAGAAGATTAGCATGGGCCCTGCGC IIIIIIIIIIIIIIIIII9IIIIIIIIII9I NH:i:1 HI:i:1 AS:i:57 nM:i:0 CB:Z:TTCTGGCGCAACTATGCAAT UR:Z:TGACTTAT UY:Z:JJJJJJJJJJ UB:Z:TGACTTAT CGGAAGACTTTATCATGATC,AATCGGTC,LH00403:55:223CL5LT1:1:1122:8193:22630 163 1 3179227 255 8S41M = 3179227 41 GACTTACATATGTTAAACCATCTCTGCATCCTTGAAATAAAACACACTT IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII NH:i:1 HI:i:1 AS:i:65 nM:i:8 CB:Z:CGGAAGACTTTATCATGATC UR:Z:AATCGGTC UY:Z:JJJJJJJJJJ UB:Z:AATCGGTC CGGAAGACTTTATCATGATC,AATCGGTC,LH00403:55:223CL5LT1:1:1122:8193:22630 83 1 3179227 255 8S41M = 3179227 -41 GACTTACATATGTTAAACCATCTCTGCATCCTTGAAATAAAACACACTT IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII NH:i:1 HI:i:1 AS:i:65 nM:i:8 CB:Z:CGGAAGACTTTATCATGATC UR:Z:AATCGGTC UY:Z:JJJJJJJJJJ UB:Z:AATCGGTC ATTGAGGATACCGATTCGCA,GCATTTAC,LH00403:55:223CL5LT1:2:2160:30282:22646 99 1 3277256 255 9S19M = 3277256 19 CGGTGAGTCTCTTAATGGAAAAAATAAA IIIIIIIIIIIIIIIIIIIIIIIIIIII NH:i:1 HI:i:1 AS:i:37 nM:i:0 CB:Z:ATTGAGGATACCGATTCGCA UR:Z:GCATTTAC UY:Z:JJJJJJJJJJ UB:Z:GCATTTAC ATTGAGGATACCGATTCGCA,GCATTTAC,LH00403:55:223CL5LT1:2:2160:30282:22646 147 1 3277256 255 9S19M = 3277256 -19 CGGTGAGTCTCTTAATGGAAAAAATAAA IIIIIIIIIIIIIIIIIIIIIIIIIIII NH:i:1 HI:i:1 AS:i:37 nM:i:0 CB:Z:ATTGAGGATACCGATTCGCA UR:Z:GCATTTAC UY:Z:JJJJJJJJJJ UB:Z:GCATTTAC CAGCCATGATTTACCGAGGC,GGTGTTTG,LH00403:55:223CL5LT1:1:1216:31253:3532 83 1 3289229 255 18M8S = 3155985 -133262 AAAAGACAAAAACCCTAAAAAACTCC 9I-9-9-I9-9--I--I-III-9--- NH:i:1 HI:i:1 AS:i:139 nM:i:5 CB:Z:CAGCCATGATTTACCGAGGC UR:Z:GGTGTTTG UY:Z:JJJJJJJJJJ UB:Z:GGTGTTTG

and the first few lines of my barcode:

AAGGCGTCGCAATCCATCTT AAGGCGTCGCACTTAACCTT AAGGCGTCGCAGGTTAGTTC AAGGCGTCGCATCATATTAG AAGGCGTCGCATCTCCGAAC AAGGCGTCGCCAACCGCTAA

And also my generated GRCm39_112 file for reference, just in case something's wrong there. GRCm39_112.csv

Any insight would be welcome.

ThuyTien1 commented 4 months ago

@msaland Can you share content of the output file/barcode_index.csv?

msaland commented 4 months ago

Hello, sorry for the late reply, here's what the top lines of the index look like: barcode_index.csv

ThuyTien1 commented 4 months ago

@msaland The files look fine to me. May I ask how large your bam file is? Can you check if there is any read in BAM file has gene information?

import pysam
bam = pysam.AlignmentFile(bamfile, "rb")
for r in  bam.fetch(until_eof=True):
    if line.has_tag("GX"):
       print(line)
       break
msaland commented 4 months ago

@ThuyTien1 ~1.6 Gb; none of the reads have gene information (it's not giving me any output for your script).

I had run the pipeline for a different dataset and that one worked (it doesn't have GX/gene information tags either).

ThuyTien1 commented 4 months ago

@msaland It's impossible to get any result if there is no GX information in any read. It's a bit weird to me that no read is confidently assigned to any gene. Is there any problem with the quality or reference annotation or sequences are not in any annotated gene?

msaland commented 4 months ago

Oh, now that you mention it, I think that was the issue. The other dataset I was working with did have GX barcodes for some of the reads.

For anyone else who runs into this specific issue: I was running STAR on the dataset that failed; STAR does not add the GX barcodes. You need to run STARsolo if you want to have the GX barcodes needed.