3
$\begingroup$

I am currently doing an analysis requiring a multiple-input script. I have a directory filled with 900+ samples.

Each sample is comprised of two-files in a fastq format. This is an example of a sample:

MIP_REV_BC_01_01A_S1_R1_001.fastq.gz 
MIP_REV_BC_01_01A_S1_R2_001.fastq.gz

The script has a Read1 & Read2 variable, and I need to make it so the R1 file is inputted for Read1 and the R2 file is inputted for Read2 for every sample. To be clear, the script is supposed to be executed on every sample separately.

This was an attempt:

with open(filename) as file:
    lines=file.readlines()
    for line in lines:
        parts[]=line.split("_")
        if parts[7]="R1":
            r1[].append(line)
        else if parts[7]="R2":
             r2[].append(line)

We were bogged down by syntax errors, but either way, this script doesn't seem to do what we wanted it to. However, it captures the idea of what we're trying to achieve.

Question I want to input a directory into the function so it automatically runs the script on the samples in the directory.

Any idea what I should do?

$\endgroup$
1
  • $\begingroup$ is the idea that you want to process both files at the same time and not store them in memory? $\endgroup$ Mar 31, 2022 at 4:47

1 Answer 1

2
$\begingroup$

Edited with explanation below. I don't really understand the objective, because I initially thought you were editing the content of fastaq file, but you seem to be editing a file which contains fastaq filenames. Is that correct? I know the syntax but I don't know really what the goal is.

Note I'm using the pathlib which has a glob function, do consider rglob under pathlib.

import re
import pathlib as Path

# pathlib is super, good to habit to use it
mypath = Path('/Users/fewbugs/workingd')
allfilenames = mypath.glob("*.fastaq.gz")
r1, r2 = [], []
x = re.compile('_')
for filename in allfilenames:
    # the 'r' is super important 
    with open(filename, 'r') as file:
    # there is no need for a readlines() 
        for line in file:
    # split via re, it is better to compile (faster) 
        parts=x.split(line)
    # you must use '=='. REMEMBER lists start at 0 not 1
           if parts[6]=="R1":
    # if you want to append to a list, declare the list first 
            r1.append(line)
    # there is no such thing as 'else if' REMEMBER lists start at 0 not 1
            elif parts[6]=="R2":
                 r2.append(line)
    # when testing good idea to trap the stuff you're missing 
            else:
                 print ("This one escaped %s" % (line))

Please be CAREFUL with my spacing, everything needs to be multiples of 4. I've not checked the code.


Forget this bit I thought you were processing the content of the fastaq file.

outfile1 = Path(mypath.parent, 'outfile1')
outfile2 = Path(mypath.parent, 'outfile2')

with open (outfile1, 'w') as f:
    for item1 in r1:
        f.write("%s\n" % item1)
with open (outfile2, 'w') as f2:
    for item2 in r2:
        f2.write("%s\n" % item2)
 

A simpler way?

... it looks like you're just processing filenames...

Might I suggest just globing no need for all this stuff.

mypath = Path('/Users/fewbugs/workingd')
r1filenames = mypath.glob("*_R1_*.fastaq.gz")
r2filenames = mypath.glob("*_R2_*.fastaq.gz")

If that is the aim, thats it! Just be careful about the glob versus rglob.

Idea If you want two paired lists this might be better done in a nested comprehension, e.g. a counter could be inserted at the central glob (*).

That may very well be what you're after. i.e.

R1 = ['stuff_R1_001.fastaq.gz', 'stuff_R1_002.fastaq.gz' 'stuff_R1_003.fastaq.gz']
R2 = ['stuff_R2_001.fastaq.gz', 'stuff_R2_002.fastaq.gz' 'stuff_R2_003.fastaq.gz']
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.