Interpreting outputs

[MadeleineOman]

Hello,

Currently we are troubleshooting the pipeline and to make sure we feel confident in our calls (that our library prep worked) we'd like to be able to visualize the mutations and read bundles better. Is there any way to get which read bundle a given mutation was called from?

Thank you !

[fa8sanger]

Hi, I am not sure I understand. In the final VCF it should indicate which read bundle supported which mutation.

[MadeleineOman]

Ah I see, the results.muts.vcf.gz file? That file takes this format:

##fileformat=VCFv4.2
##source=NanoSeq pipeline
##FILTER=<ID=PASS,Description="All filters passed">
##reference=file://raw_data/sequencing/reference/genome.fa
##contig=<ID=JF411744.1,length=330611>
##ALT=<ID=*,Description="Represents allele(s) other than observed.">
##INFO=<ID=TRI,Number=1,Type=String,Description="Pyrimidine context, trinucleotide substitution">
##INFO=<ID=BBEG,Number=1,Type=String,Description="Read bundle left breakpoint">
##INFO=<ID=BEND,Number=1,Type=String,Description="Read bundle right breakpoint">
##INFO=<ID=QPOS,Number=1,Type=Integer,Description="Read position closest to 5-prime end">
##INFO=<ID=DEPTH_FWD,Number=1,Type=Integer,Description="Read bundle forward reads depth">
##INFO=<ID=DEPTH_REV,Number=1,Type=Integer,Description="Read bundle reverse reads depth">
##INFO=<ID=DEPTH_NORM_FWD,Number=1,Type=Integer,Description="Matched normal forward reads depth">
##INFO=<ID=DEPTH_NORM_REV,Number=1,Type=Integer,Description="Matched normal reverse reads depth">
##FILTER=<ID=dbsnp,Description="Common SNP site">
##FILTER=<ID=shearwater,Description="Noisy site">
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
JF411744.1	35142	.	C	T	.	PASS	TRI=ACA>T;BBEG=35052;BEND=35161;QPOS=19;DEPTH_FWD=6;DEPTH_REV=4;DEPTH_NORM_FWD=2002;DEPTH_NORM_REV=1554
JF411744.1	35174	.	G	A	.	PASS	TRI=ACC>T;BBEG=34988;BEND=35258;QPOS=84;DEPTH_FWD=5;DEPTH_REV=2;DEPTH_NORM_FWD=1575;DEPTH_NORM_REV=1695
JF411744.1	35293	.	G	A	.	PASS	TRI=CCA>T;BBEG=35246;BEND=35452;QPOS=47;DEPTH_FWD=5;DEPTH_REV=3;DEPTH_NORM_FWD=1827;DEPTH_NORM_REV=1763

I'm assuming the "Read bundle left breakpoint" and "Read bundle right breakpoint" are the "start" and "end" respectively in the RB id tag (format RB:Z:chrom,start,end,left_tag,right_tag ). Is this the case? If so, is there no way to differentiate b/t read bundles that have the same coordinates but different tags? (Happens occasionally in my dataset).

Another question, how can I get all the read bundles that were identified from the pipeline? I see read bundles are listed in the Sample.RBs file output by the efficiency_nanoseq.pl script (attached), however the ones listed correspond to the "total # RBs" as outlines in the sample.tsv summary output by the efficiency script (attached). How can I get the total number of "OK" RBs? I tried filtering for lines with >=2 in both columns ( I assume the forward and reverse strand depth in the .RBs file) but that didn't quite align with the sample.tsv summary.

Thank you!

DL001.tsv

DL001_header.RBs

[fa8sanger]

Hi Madeleine, it seems you are using an older version of the pipeline. The current one also outputs a field named "BTAG" that contains the other elements required to reconstruct the whole read bundle id. You may need to run the last step of the pipeline again with a more recent version
These are the current tags:
##INFO=<ID=BTAG,Number=1,Type=String,Description="Read bundle tag (duplex barcode)">
##INFO=<ID=BBEG,Number=1,Type=String,Description="Read bundle left breakpoint">
##INFO=<ID=BEND,Number=1,Type=String,Description="Read bundle right breakpoint">
##INFO=<ID=TRI,Number=1,Type=String,Description="Pyrimidine context, trinucleotide substitution">
##INFO=<ID=QPOS,Number=1,Type=Integer,Description="Read position closest to 5-prime end">
##INFO=<ID=DEPTH_FWD,Number=1,Type=Integer,Description="Read bundle forward reads depth">
##INFO=<ID=DEPTH_REV,Number=1,Type=Integer,Description="Read bundle reverse reads depth">
##INFO=<ID=DEPTH_NORM_FWD,Number=1,Type=Integer,Description="Matched normal forward reads depth">
##INFO=<ID=DEPTH_NORM_REV,Number=1,Type=Integer,Description="Matched normal reverse reads depth">
##INFO=<ID=DPLX_ASXS,Number=1,Type=Integer,Description="AS-XS for duplex">
##INFO=<ID=DPLX_CLIP,Number=1,Type=Integer,Description="Clipping for duplex">
##INFO=<ID=DPLX_NM,Number=1,Type=Integer,Description="Mismatches in duplex">
##INFO=<ID=BULK_ASXS,Number=1,Type=Integer,Description="AS-XS for bulk">
##INFO=<ID=BULK_NM,Number=1,Type=Integer,Description="Mismatches in bulk">

[fa8sanger]

About your other question, those results are only for JF411744.1. If you want the pipeline to output all of them for the entire genome, you would need to modify the efficiency script.
Regarding the discrepancy that you report between the xxx.RBs file and the output in the efficiency results, I am not sure, they should match. Indeed, I had a look at some examples and they match. Ah, I see, you have $2 < 1 rather than $2 > 1 :-)