Removing reads with "N" > 4

[marcela]

Hi!

I was wondering if there is a tool to remove reads that have a specific number of "N's". Fastx has the fastx_trimmer, but it takes positions to trim, without looking at the nucleotide content. Any suggestions?

Thanks!

[NicoBxl]

in R it's easy

countN <- function(a){
  return(sum(strsplit(a,"")[[1]]=="N")==4)
}

a <- "ATGCNNNN"
b <- "ATGCNN"
x <- c(a,b)
sapply(x,countN)

ATGCNN   ATGN 
     TRUE      FALSE

[marcela]

Quite simple! Thanks!

[NicoBxl]

you're welcome

[schmima]

hm - don't know of any - but is relatively easy task - heres a quick and dirty script (assumes that your infile has no header) - most probably not too fast due to simultanous reading and writing. But it should work (did not test - was originally something else):

[EDIT - THE TABS ARE NOT CORRECTLY PRINTED HERE... DONT KNOW IF THEY WILL APPEAR IF YOU COPY IT]

############

import sys

if len(sys.argv) < 5:
print 'usage:\n\tpython %s infile outfile char maxnum' % (str(sys.argv[0]))
sys.exit(0)

in_file = sys.argv[1]
out_file = sys.argv[2]
char = sys.argv[3]
maxnum = int(sys.argv[4])

infile = open(in_file, 'r')
outfile = open(out_file, 'w')

prevLines = []
for line in infile:
prevLines.append(str(line[:-1]))
if len(prevLines) == 4:
nucID = prevLines[0]
nucSeq = prevLines[1]
qualID = prevLines[2]
qualSeq = prevLines[3]
if nucSeq.count(char) <= maxnum:
outstring = '\n'.join([nucID,nucSeq,qualID,qualSeq])+'\n'
outfile.write(outstring)
prevLines = []

if len(prevLines) == 4:
charcounter = 0
nucID = prevLines[0]
nucSeq = prevLines[1]
qualID = prevLines[2]
qualSeq = prevLines[3]
if nucSeq.count(char) <= maxnum:
outstring = '\n'.join([nucID,nucSeq,qualID,qualSeq])+'\n'
outfile.write(outstring)


outfile.close()
infile.close()


##############

copy the part between the ###### into an empty textfile, save as XY.py and run:
python XY.py pathtoinfile pathtooutfile characteryouwanttoremove threshold

eg
python XY.py /home/me/myreads.fastq /home/me/myreads_filtered.fastq N 4

(removes all reads with more than 4 Ns - ie writes the others to outfile)

[schmima]

ad tabs: if necessary use this code here and replace every * with one tab. ** are double tabs and so on:

import sys

if len(sys.argv) < 5:
*print 'usage:\n\tpython %s infile outfile char maxnum' % (str(sys.argv[0]))
*sys.exit(0)

in_file = sys.argv[1]
out_file = sys.argv[2]
char = sys.argv[3]
maxnum = int(sys.argv[4])

infile = open(in_file, 'r')
outfile = open(out_file, 'w')

prevLines = []
for line in infile:
*prevLines.append(str(line[:-1]))
*if len(prevLines) == 4:
**nucID = prevLines[0]
**nucSeq = prevLines[1]
**qualID = prevLines[2]
**qualSeq = prevLines[3]
**if nucSeq.count(char) <= maxnum:
***outstring = '\n'.join([nucID,nucSeq,qualID,qualSeq])+'\n'
***outfile.write(outstring)
**prevLines = []

if len(prevLines) == 4:
*charcounter = 0
*nucID = prevLines[0]
*nucSeq = prevLines[1]
*qualID = prevLines[2]
*qualSeq = prevLines[3]
*if nucSeq.count(char) <= maxnum:
**outstring = '\n'.join([nucID,nucSeq,qualID,qualSeq])+'\n'
**outfile.write(outstring)


outfile.close()
infile.close()