lh3 / psmc

Implementation of the Pairwise Sequentially Markovian Coalescent (PSMC) model
Other
146 stars 60 forks source link

Empty .psmcfa file #32

Open JuliaLopezDelgado opened 3 years ago

JuliaLopezDelgado commented 3 years ago

I am running the following code:

samtools mpileup -uf ref.fasta sorted.bam | bcftools call -c | vcfutils.pl vcf2fq | gzip > Cmar.fq.gz
psmc/utils/fq2psmcfa -q20 Cmar.fq.gz > Cmar.psmcfa 
psmc/psmc -N25 -t15 -r5 -p "4+25*2+4+6" -o Cmar.psmc Cmar.psmcfa

I don't get any errors when running it but it produces an empty .psmcfa file, and consequently a .psmc file with NAN values.

Could you please help me with this issue?

Thanks in advance, Julia

JuliaLopezDelgado commented 3 years ago

I have tried running fq2psmcfa with different parameters (i.e. lower q values, and small window and scaffold sizes) but I still get an empty .psmcfa file.

plubbe commented 3 years ago

Try checking the contents of each intermediate file- your input . bam and the .fq.gz output. I had the same problem which was solved when I realized I was pointing to empty files due to a typo in my paths. Hopefully it will be an easy solution like that.

hennelly commented 3 years ago

Hi All,

I'm also getting the same issue, where my script works for the conversion of the bam to fq.gz files, however it ends up producing empty .psmcfa files. There's no information in the error file, so its been difficult to trouble shoot. Hope someone can help with this issue as well.

samtools mpileup -q 30 -Q 30 -C50 -uf ref.fa sample_sorted_proper_nodups_readgroup.bam | bcftools call -c | vcfutils.pl vcf2fq -d 5 -D 30 | gzip > sample.fq.gz

/bin/psmc/utils/fq2psmcfa -q30 sample.fq.gz > sample.psmcfa

Thanks,

-Lauren

gubrins commented 3 years ago

Hi all,

I was having the same problem and I discovered that my bam file was corrupted, so as indicated before, it could be related to some kind of problems with your input files.

Biscuite-wzy commented 2 years ago

I was having the same problem

hennelly commented 2 years ago

Hi Everyone,

I think I figured out my issue. I used fastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) to check the quality of my fastq files, which were downloaded from NCBI SRA.

I then found that the minimum mapping quality of the samples that were not working for PSMC were on average ~ 25, which in my script I was filtering for 30 using the -Q flag in samtools mpileup. I changed the filter to -Q 20, and the scripts seem be working.

linhai1221 commented 2 years ago

if ((double)n_good_bases / len >= 0.2 && n_good_bases >= n_min_good) change 0.2 to 0

Sixth-Grandpa commented 2 years ago

if ((double)n_good_bases / len >= 0.2 && n_good_bases >= n_min_good) change 0.2 to 0

Hi, linhai1221. I have tried this change. But the psmcfa file was still empty. I look forward to hearing from you.

Sixth-Grandpa commented 2 years ago

Hi Everyone,

I think I figured out my issue. I used fastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) to check the quality of my fastq files, which were downloaded from NCBI SRA.

I then found that the minimum mapping quality of the samples that were not working for PSMC were on average ~ 25, which in my script I was filtering for 30 using the -Q flag in samtools mpileup. I changed the filter to -Q 20, and the scripts seem be working.

Thanks for your suggestion. As your method, I got the correct psmcfa file!

Candce commented 1 year ago

Hi I cannot seem to get a psmcfa file for the example file using the code utils/fq2psmcfa -q20 diploid.fq.gz > diploid.psmcfa. It seems as if psmc isnt working and I was wondering if its happened before. I went to utils and said make but it doesnt look like its compiling properly. I have tried with my own data and cannot get a psmcfa file using the following code utils/fq2psmcfa /n/holyscratch01/edwards_lab/Lachenicht/bamOUTsex/r.fq > r.psmcfa. I also found the same tools pile up worked as the file isnt empty. I have tried changing the q flag but it looks like the problem is with the software. Did I miss something in the compile step? thanks Candice