For loop for aligning multiple paired end fastq files
1
0
Entering edit mode
9 weeks ago
kcarery • 0

Hey all,

Could you advise what I am doing wrong?

Files looks like this:

OV63.FCD1U51ACXX_L7_IACTTGA.fastq_R1.fastq
OV63.FCD1U51ACXX_L7_IACTTGA.fastq_R2.fastq


I tried this code to align to reference genome:

#Aligning to reference genome
for f in *.fastq | rev | cut -c 10- | rev | uniq
do
STAR --genomeDir ~/file/path/Indexes/ncbi-genomes-2022-09-19/
--readFilesIn ${f}_R1.fastq {f}_R2.fastq --outSAMtype BAM Unsorted --outReadsUnmapped Fastx --outFileNamePrefix AlignedHG38/${f}
done


I get these errors:

-bash: syntax error near unexpected token |'
Nov 29 13:25:06 ...... FATAL ERROR, exiting

star • 573 views
0
Entering edit mode
9 weeks ago
Ram 37k

The shell stops processing your pipe at *.fastq. If you need to use something like the old version of my bash loop command, you need to use a subshell to do the piped command: $(ls *.fastq | rev | cut -c 10- | rev | uniq) instead of simply *.fastq | rev | cut -c 10- | rev | uniq. Instead of using the hardcoded cut -c 10- (which will definitely fail in your use case), adapt the updated $(ls *R1*.fastq | sed -r 's/_R1_001[.]fastq//' | uniq) - your FASTQ files don't have the _001 part. Again, for details, see my older post.

0
Entering edit mode

Thanks for your response! I read your old version link and I understand why sed is better. I still am confused how do you add in the R1 and R2, I am trying to align it to a reference but this would only take in to account R1 right? I am using STAR aligner. for the STAR --readfilesin I must list both R1 and R2

How would I also add in R2?

#Aligning to reference genome

for f in $(ls *R1*.fastq | sed -r 's/_R1[.]fastq//' | uniq) do STAR --genomeDir ~/file/path/Indexes/ncbi-genomes-2022-09-19/ --runThreadN 20 --readFilesIn${f}_R1.fastq {f}_R2.fastq --outSAMtype BAM Unsorted --outReadsUnmapped Fastx
done

0
Entering edit mode

I think there's either another gap in understanding or a typo: you're using the R1 file variable right (${f}_R1.fastq) but are missing the $ sign in the R2 (it should be ${f}_R2.fastq, not {f}_R2.fastq). Do you understand (a) what the $ does, and (b) what the {} after the $ are for? ADD REPLY 0 Entering edit mode Yes I see I was missing that. I didn't even catch it. What hours at the computer will do to the mind... The $ sign tells us to call that variable correct? Similar to 3x+5 type of thing? Why the function of x is....in this case would be this $(ls *R1*.fastq | sed -r 's/_R1[.]fastq//' | uniq) ....that's how I understood it... ADD REPLY 0 Entering edit mode $ sign tells us to call that variable

Close. $ tells the interpreter to look for a variable whose name is the sequence of characters following the $. So if your variable were to be named my_var (example: my_var=ABCDE_L001, then the shell would replace $my_var with ABCDE_L001 wherever it encounters $my_var.

However, if you were to use something like $my_var_R1, the shell would replace that with an empty string, as it doesn't know a variable named my_var_R1 and it does not know that you intend to append _R1 to $my_var. This is where the {} come into play. They "isolate" the variable name from the rest of the string, and using them is a good practice in general. For example, the shell would understand {my_var}_R1 and you'd see ABCDR_L001_R1. For even more info, see my comment from that thread: bash loop for alignment RNA-seq data ADD REPLY 0 Entering edit mode Ahhh I see. Okay, that makes sense. Its the prefix.... okay cool. So I did refer to the thread but still not sure what I am missing. I attempted to also add in a echo command of the STAR alignment because the issue is no longer about being able to obtain the correct file. I was able to do that by practicing with a few files. For example: #check how it will run on files with the full code for f in(ls *R1*.fastq | sed -r 's/_R1[.]fastq//' | uniq)
do
echo STAR --genomeDir ~/file/path/Indexes/ncbi-genomes-2022-09-19/
echo --readFilesIn ${f}_R1.fastq${f}_R2.fastq
echo --quantMode TranscriptomeSAM GeneCounts
echo --twopassMode Basic
echo --chimOutType Junctions
echo --outSAMtype BAM Unsorted
echo --outFileNamePrefix ../starAlignedHG38/\${f}
echo --outSAMtype BAM Unsorted
done


It gave me this:

STAR --genomeDir /file/path/Indexes/ncbi-genomes-2022-09-19/
--quantMode TranscriptomeSAM GeneCounts
--twopassMode Basic
--chimOutType Junctions
--outSAMtype BAM Unsorted
--outFileNamePrefix ../starAlignedHG38/OV88.FCC1PU8ACXX_L3_IGGCTAC.fastq_
--outSAMtype BAM Unsorted

STAR --genomeDir /file/path/Indexes/ncbi-genomes-2022-09-19/
--quantMode TranscriptomeSAM GeneCounts
--twopassMode Basic
--chimOutType Junctions
--outSAMtype BAM Unsorted
--outFileNamePrefix ../starAlignedHG38/OV89.fastq_
--outSAMtype BAM Unsorted


Which is what I would want to see.... It is still giving me the error of it can read the ReadFilesIn tab, saying it can not read it and I am confused.... I've followed you all way of doing it and checked it with echo

0
Entering edit mode

Please use the formatting bar (especially the code option) to present your post better. You can use backticks for inline code (text becomes text`), or select a chunk of text and use the highlighted button to format it as a code block. If your code has long lines with a single command, break those lines into multiple lines with proper escape sequences so they're easier to read and still run when copy-pasted. I've done it for you this time.

0
Entering edit mode

That is odd. Can you copy-paste the exact error message again, please? The smallest of differences in the error message can help us crack the problem.

EDIT:

Try running the command without breaking it into multiple lines. You're not using escape sequences and STAR is probably reading only the first parameter, stopping at the newline after it.

0
Entering edit mode

Thank you for your support, it is finally working

0
Entering edit mode

What did you have to do to get it working?