I am introducing SNVs into specific samples in order to estimate false negative rates for a variant calling pipeline. I know reads can be simulated but I would actually prefer to use the real data so as to keep everything else equal.
Given a read ID, I want to edit a single basecall (e.g. the 12th base) for just that read from within a large FASTQ file containing millions of reads.
example, i want to change the 12th base ('C') in read 31027 to a 'T':
@70630 1:N:0:ATCACG GAAGGTCCATGGATAATACTCAATTTTCCACAACAGCTTTTGTACTCTAGATCATTGATATTTACCAAAAGTCACTCAAACTCATCCTATGCATAATTCTAGTCCACCAATCATGATATGATGGAGAACATGGTTGTAATCAGGAAGACAG + DDDDDHIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIHIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIGHHIIIIIIIIIIIIIIIIIIIIIIHIIHIHIIIHHIGIIIIHHHHIIH@ @31027 1:N:0:ATCACG CAAAAGTCACTCAAACTCATCCTATGCATAATTCTAGTCCACCAATCATGATATGATGGAGAACATGGTTGTAATCAGGAAGACAGATAAAGCAGCAGACCAAAAGTAATCTGAGAAATTATATTTGAATCACTCAGATATACATCAAATA + DDBDDHC@GHFHHIIE@CEHHIHEFFCGFCH?HHEGFCCEHHCGH@HFHHHCHHIIIHEHIIHII@CECHIHIIIEECCDGEHFFHHHHEHHFHHGGIHHHDGHFHIIIGHHHHHHEHHIIHIIIGGHHHCHHHCHIIIHHIEH@GHIHIC @87319 1:N:0:ATCACG CAATTAAGCTTTGGCAACGGTGGTCAAGATGAGATGCATATGGAGATAATAACTAAAAGTCAATCGAGACTCATCGTATGCATATTTCTAGTCCATCGATCATGAAATGATAGGATAGCTAGAATGAAAAGTAAATTTCCAGAAGGTCCAT + DDDDDIIIIIIIIIIIIIHIHHIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIHHIIIIIIIIIIIIIIIIIHIIHIIIIIIIIIIIIHIIIIIIIIIHIHIH
Of course I could stream through the entire FASTQ file rewriting everything to a new file until i get to the read of interest, do the edits to that read, write it out to the new file, and then keep going until the end. BioPython would work fine for this as would other approaches.
However, is there an efficient way to do this without the full read/write stream? Could this be done with sed/awk ? What about an indexed fastq file?
cheers David