Issue with reverting bam file back to fastq files

Question

[williamsbrian5064]

Hi,

I have some old bam files that I don't have the original fastq files for. I tried using picard Samtofastq but I am having an issues when I realign the reads. I have been using picard to revert the old bam files back to fastq files, bwa to align the reads, and picard to sort and deduplicate. I used samtools flagstat to check out what was going on with my new bam files and this is what I got

473763117 + 0 in total (QC-passed reads + QC-failed reads)
1618429 + 0 secondary
0 + 0 supplementary
24690025 + 0 duplicates
473247331 + 0 mapped (99.89% : N/A)
472144688 + 0 paired in sequencing
236072344 + 0 read1
236072344 + 0 read2
466910892 + 0 properly paired (98.89% : N/A)
471430496 + 0 with itself and mate mapped
198406 + 0 singletons (0.04% : N/A)
3805718 + 0 with mate mapped to a different chr
2458472 + 0 with mate mapped to a different chr (mapQ>=5)

I have checked out my files and it doesn't look like it is an issue with bwa or picard sort/deduplicate. The only thing I can think of is that there is some issue with picard samtofastq

This is the command I use

java -jar picard.jar SamToFastq I=C01767.bam FASTQ=C01767.R1.fastq SECOND_END_FASTQ=C01767.R2.fastq

Am I doing something wrong here?

I know GATK has a way to do this, but their method only give you one files here. I kind of want both of the fastq files

[Eric Lim]

Never used picard to convert bam to fastq. For paired end, I use bedtools. Did you sort the bam by name before conversion?

[williamsbrian5064]

I did not sort the bam files by name before conversion. What does that do?

[GenoMax]

That brings paired-end reads (R1/R2) next to each other in the alignment file.

[GenoMax]

Why not use samtools fastq? (equivalent to samtools bam2fq).

[williamsbrian5064]

I can give it a try. Hopefully that will help

[h.mon]

What is the problem? Your alignment stats seems fine to me.

[WouterDeCoster]

Maybe I missed it, but what's wrong?

[williamsbrian5064]

When I use samtools mpileup, there are no variants being called. All the other samples have variants being called except the samples that were originally old bam files.

[Devon Ryan]

What is the exact command you're using. samtools mpileup doesn't directly call variants any more, bcftools call does.

[WouterDeCoster]

I think that this is important information you should mention in your initial question, rather than to make us guess what the question is actually about.

[finswimmer]

Hello williamsbrian5064,

you needn't reply to each comment with the same text. Everyone who is participating in this thread gets noticed if there is a new post. So we decided to delete your duplicate posts to keep the focus on your problem.

fin swimmer

[williamsbrian5064]

Thanks for doing that finswimmer! I apologize for that. I didn't realize that everyone got a notification.

Sorry that I didn't mention my initial question. This is the command that I am putting in. Some of the tags are a bit out of data. I mean samtools mpileup works fine with the other 30 samples, so I don't know if it is a samtools mpileup vs bcftools call issue?

samtools mpileup --redo-BAQ --min-BQ 30 --per-sample-mF -C 50 --output-tags DP,DV,DPR,INFO/DPR,DP4,SP -r X:1-126427096 -u -f inputfasta.fa --BCF input.bam > output.vcf

[GenoMax]

It would be helpful to edit the original question and add this information there. This thread has already got somewhat out of hand.

[finswimmer]

Your flagstats looks fine to ne. What do you think is wrong?

fin swimmer