Redirect alignments to file

[Pkaps25]

Would it be possible to redirect the alignment output of danbing-tk align to a file, perhaps through an argument to the align command? Right now, the only way I see to get reads is to grep from stdout, which becomes difficult when running jobs in parallel using snakemake.

[joyeuxnoel8]

I usually do this by changing the command in snakemake to something like danbing-tk -a <other_arguments> | gzip > out.aln.gz.

[joyeuxnoel8]

Ah good question, the locus (0-based) denotes the order in pan.tr.mbe.v2.bed.

[Pkaps25]

Thank you. If I were to make a list of all loci with nonzero length estimates in a given sample and a list of loci that had an aligned read during the alignment step, should I expect these two lists to be the same? I am thinking yes, because of the length prediction formula. If the kmer count for a locus is nonzero, which occurs when a read aligns to that locus, then the length estimate should be nonzero.

[joyeuxnoel8]

Theoretically yes! But on real data I would expect noise. So measuring the correlation between the length and kmer dosage might be more appropriate.

[Pkaps25]

Thank you, I am now more interested in identifying what VNTRs are represented in my data than estimating their lengths and was hoping to skip the length prediction steps as I can't get accurate LSBs. To get LSBs, I have been using a single reference genome as a ground truth, which naturally leads to flawed LSBs. I can either run length estimation using these flawed LSBs and use loci with length>0, run length estimation using only kmer dosage and use loci with length >0, or use loci that were aligned to in the align step. I'm unsure which option is best, do you have any suggestions?

Also, following up on your last point, why would the correlation between length and kmer count help with reducing the impact of noise? When using kmer dosage, one of my samples has length estimated for ~1500 loci/15000, but when using a count of unique loci from the alignment output, I am getting ~5000 unique loci. I wonder if this has something to do with the point you made.

[joyeuxnoel8]

I'm unsure which option is best, do you have any suggestions?

I might need to know your overall goal better to make suggestions. Are you trying to compare the differences between TB strains by examining their VNTR contents? Let me know.

Also, following up on your last point, why would the correlation between length and kmer count help with reducing the impact of noise?

I was assuming that your goal is to count the number of loci with nonzero kmer counts (let it be N1). In this scenario, N1 would be sensitive to any misalignments between loci. I'm not very familiar with TB VNTRs, but if there exists loci that are highly similar to each other, then misalignments are inevitable.

When using kmer dosage, one of my samples has length estimated for ~1500 loci/15000, but when using a count of unique loci from the alignment output, I am getting ~5000 unique loci. I wonder if this has something to do with the point you made.

The alignment output also contains mapping to flanking regions. If you're only interested in VNTR kmer counts, you should stick to the tr.kmers output.