Hello,

I am looking to read every sequence of a fastq file, insert a specific set of nucleotides at the beginning of each sequence. Then save that file into a new fastq file.

I am using the shortRead package for that, but am not sure how to overwrite the new sequence. This is what I have so far.

library("ShortRead")
input_1 <- readFastq("Sample_1.fastq") #input forward read

output_1 <- input_1 #making a copy 

insertNuc <- "TTTTT"  #to be inserted into sequence
iCount = 1
while(iCount <= length(input_1)) {
    seq1 <- sread(input_1)[iCount]  #read sequence one by one
    seq1New <- paste(insertNuc,seq1,sep="")  #insert nucleotides
    sread(output_1)[iCount] <- seq1New  # write back new sequence to fastq file  ????

    iCount <- iCount + 1
}
# write fastq
writeFastq(output_1, destination, "output_1") 

I need help with

1)  write back new sequence to fastq file in same order (basically I need to replace old sequence with new sequence)

2)  I am using loops to read through entire fastq file. If anyone knows a faster way to do this, with apply() or anything other packages , I would appreciate any tips on this.

Thanks, K

Don't read the records one-by-one, but in 'chunks' of say 1000000. Do this using ?FastqStreamer; the example on this help page provides some illustration of how to do this using repeat {}; see also?GenomicFiles::reduceByYield.

For example, this sets up a stream to read 100 records at a time (as mentioned, work on much larger chunks), and then extracts the first 100 records

example(FastqStreamer)         # defines 'fl'
strm = FastqStreamer(fl, 100)  # much larger chunk, e.g,. 1000000
chunk = yield(strm)

Manipulate the components of the resulting object, maybe using ?replaceAt, and then re-create the object

sr = sread(chunk); sr = replaceAt(sr, 1, "TTTT")
qual = quality(quality(chunk)); qual = replaceAt(qual, 1, "    ")
fq = ShortReadQ(sr, do.call(class(quality(chunk)), list(qual)), id(chunk))

And write the result

writeFastq(fq, "updated.fastq", mode="a")

A complete example is

library(ShortRead)
library(GenomicFiles)
example(FastqStreamer)  # for the path 'fl' to a fastq file

YIELD = yield
MAP = function(x, fout) {
    sr = sread(x); sr = replaceAt(sr, 1, "TTTT")
    qual = quality(quality(x)); qual = replaceAt(qual, 1, "    ")
    x <- ShortReadQ(sr, do.call(class(quality(x)), list(qual)), id(x))
    writeFastq(x, fout, "a")
    length(x)
}
fout = tempfile()
reduceByYield(FastqStreamer(fl), YIELD, MAP, fout=fout)

With 

> readFastq(fl)
class: ShortReadQ
length: 256 reads; width: 36 cycles
> readFastq(fout)
class: ShortReadQ
length: 256 reads; width: 40 cycles