This is a test version of Biostars. For the public version, visit https://www.biostars.org.
shell question: How to run multiple paired end alignment in shell?

Hi there,

the following command works with bash(Mac command line), but not with unix shell(after entering screen command), how do I rewrite it?

for i in $(ls ~/mbd1rnaseq/scfastq | rev | cut -c 9- | rev | uniq)

do  rsem-calculate-expression -p 8 --forward-prob 0.95 --bowtie2 --bowtie2-path /mnt/software/bowtie2/default --output-genome-bam   --paired-end ~/mbd1rnaseq/scfastq/${i}R1.fastq ~/mbd1rnaseq/scfastq/${i}R2.fastq ~/mbd1rnaseq/mm9.rsem ${i}

done |& tee terminaloutput.txt

Basically , file input names got changed.

Error: Must specify at least one read input with -U/-1/-2
(ERR): bowtie2-align exited with value 1
sh: 32m05_CGTACTAG-GCGTAAGA_R2.fR1.fastq: command not found
sh: 32m05_CGTACTAG-GCGTAAGA_R2.fR2.fastq: command not found
sh: 32m05_CGTACTAG-GCGTAAGA_R2.f.temp/o such file or directo

my file names are like this:

11_AGGCAGAA-TATCCTCT_R1.fastq  39_CGTACTAG-ACTGCATA_R1.fastq  60_GGACTCCT-TATCCTCT_R1.fastq  77_CTCTCTAC-AGAGTAGA_R1.fastq
11_AGGCAGAA-TATCCTCT_R2.fastq  39_CGTACTAG-ACTGCATA_R2.fastq  60_GGACTCCT-TATCCTCT_R2.fastq  77_CTCTCTAC-AGAGTAGA_R2.fastq
rna-seq alignment r software error
• 719 views • link
written by John John

What's the output of ls ~/mbd1rnaseq/scfastq?

• 0 views • link
written by finswimmer finswimmer 
01_TAAGGCGA-GCGTAAGA_R1.fastq  18_GGACTCCT-CTCTCTAT_R1.fastq  48_AGGCAGAA-AAGGAGTA_R1.fastq  67_TAGGCATG-CTCTCTAT_R1.fastq
01_TAAGGCGA-GCGTAAGA_R2.fastq  18_GGACTCCT-CTCTCTAT_R2.fastq  48_AGGCAGAA-AAGGAGTA_R2.fastq  67_TAGGCATG-CTCTCTAT_R2.fastq
02_TAAGGCGA-CTCTCTAT_R1.fastq  19_GGACTCCT-TATCCTCT_R1.fastq  49_AGGCAGAA-CTAAGCCT_R1.fastq  68_TAGGCATG-TATCCTCT_R1.fastq
02_TAAGGCGA-CTCTCTAT_R2.fastq  19_GGACTCCT-TATCCTCT_R2.fastq  49_AGGCAGAA-CTAAGCCT_R2.fastq  68_TAGGCATG-TATCCTCT_R2.fastq
03_TAAGGCGA-TATCCTCT_R1.fastq  20_GGACTCCT-AGAGTAGA_R1.fastq  51_TCCTGAGC-CTCTCTAT_R1.fastq  69_TAGGCATG-AGAGTAGA_R1.fastq
03_TAAGGCGA-TATCCTCT_R2.fastq  20_GGACTCCT-AGAGTAGA_R2.fastq  51_TCCTGAGC-CTCTCTAT_R2.fastq  69_TAGGCATG-AGAGTAGA_R2.fastq
04_TAAGGCGA-AGAGTAGA_R1.fastq  21_TAGGCATG-GCGTAAGA_R1.fastq  52_TCCTGAGC-TATCCTCT_R1.fastq  70_TAGGCATG-GTAAGGAG_R1.fastq
04_TAAGGCGA-AGAGTAGA_R2.fastq  21_TAGGCATG-GCGTAAGA_R2.fastq  52_TCCTGAGC-TATCCTCT_R2.fastq  70_TAGGCATG-GTAAGGAG_R2.fastq
• 0 views • link
written by John John

This is a problem with your cut -c 9- statement, because it hard-codes the number of characters to cut. Use sed with a regex instead, that way you can be sure that the right substring is picked regardless of length.

Check out the output of ls 32m05_CGTACTAG-GCGTAAGA_R2.f* to see what's in the file name that's breaking your cut.

Your loop might improve with this:

for f in $(ls -1 ~/mbd1rnaseq/scfastq/*.fastq | sed 's/R[[:digit:]].fastq//' | uniq)
do
    rsem-calculate-expression -p 8 --forward-prob 0.95 --bowtie2 --bowtie2-path /mnt/software/bowtie2/default --output-genome-bam   --paired-end ${f}R1.fastq ${f}R2.fastq ~/mbd1rnaseq/mm9.rsem $(basename $f)
done

Note: ls /some/dir gives you just basenames of the files in the dir, but ls /some/dir/*.glob gives you full relative path to each matching filename, removing the need to specify the path to each FQ file in your command.

• 0 views • link
written by Ram Ram

It gives like this RamRS

Fail to create folder /home/BIOTECH/xzhao/mbd1rnaseq/scfastq/70_TAGGCATG-GTAAGGAG_.temp.
Fail to create folder /home/BIOTECH/xzhao/mbd1rnaseq/scfastq/71_TAGGCATG-ACTGCATA_.temp.
Fail to create folder /home/BIOTECH/xzhao/mbd1rnaseq/scfastq/72_TAGGCATG-AAGGAGTA_.temp.
Fail to create folder /home/BIOTECH/xzhao/mbd1rnaseq/scfastq/73_TAGGCATG-CTAAGCCT_.temp.
Fail to create folder /home/BIOTECH/xzhao/mbd1rnaseq/scfastq/74_CTCTCTAC-GCGTAAGA_.temp.
Fail to create folder /home/BIOTECH/xzhao/mbd1rnaseq/scfastq/75_CTCTCTAC-CTCTCTAT_.temp.
Fail to create folder /home/BIOTECH/xzhao/mbd1rnaseq/scfastq/76_CTCTCTAC-TATCCTCT_.temp.
Fail to create folder /home/BIOTECH/xzhao/mbd1rnaseq/scfastq/77_CTCTCTAC-AGAGTAGA_.temp.
Fail to create folder /home/BIOTECH/xzhao/mbd1rnaseq/scfastq/78_CTCTCTAC-GTAAGGAG_.temp.
• 0 views • link
updated by Ram Ram • written by John John

You seem to lack write permission to ~/mbd2ranseq/scfastq. You can either address that or supply a different temp directory to rsem-calculate-expression using the --temporary-folder <string> option.

• 0 views • link
written by Ram Ram

Still It produces the same error.

Error: Must specify at least one read input with -U/-1/-2
(ERR): bowtie2-align exited with value 1
sh: 32m12_AGGCAGAA-AGAGTAGA_R1.fR1.fastq: command not found
sh: 32m12_AGGCAGAA-AGAGTAGA_R1.fR2.fastq: command not found
sh: 32m12_AGGCAGAA-AGAGTAGA_R1.f.bam: command not found
"/mnt/software/bowtie2/default/bowtie2 -q --phred33 --sensitive --dpad 0 --gbar 99999999 --mp 1,1 --np 1 --score-min L,0,-0.1 -I 1 -X 1000 --no-mixed --no-discordant -p 8 -k 200 -x /home/BIOTECH/xzhao/mbd1rnaseq/mm9.rsem -1 /home/BIOTECH/xzhao/mbd1rnaseq/scfastq/12_AGGCAGAA-AGAGTAGA_R1.fR1.fastq -2 /home/BIOTECH/xzhao/mbd1rnaseq/scfastq/12_AGGCAGAA-AGAGTAGA_R1.fR2.fastq | samtools view -S -b -o /home/BIOTECH/xzhao/ks/12_AGGCAGAA-AGAGTAGA_R1.f.
• 0 views • link
written by John John

I have a feeling that your ls is customized with shell colors and that could be interfering with operations, as there is no file that begins 32m* in that dir. Try using /bin/ls instead of just ls.

• 0 views • link
written by Ram Ram 
$ ls

sample1_R1.fastq  sample1_R2.fastq  sample2_R1.fastq  sample2_R2.fastq

script:

$ for i in *R1*.fastq; do echo rsem-calculate-expression -p 8 --forward-prob 0.95 --bowtie2 --bowtie2-path /mnt/software/bowtie2/default --output-genome-bam   --paired-end ~/mbd1rnaseq/scfastq/$i ~/mbd1rnaseq/scfastq/${i/%R1.fastq/R2.fastq} ~/mbd1rnaseq/mm9.rsem ${i} ${i/%R1.fastq/R2.fastq}; done


rsem-calculate-expression -p 8 --forward-prob 0.95 --bowtie2 --bowtie2-path /mnt/software/bowtie2/default --output-genome-bam --paired-end /home/john_doe/mbd1rnaseq/scfastq/sample1_R1.fastq /home/john_doe/mbd1rnaseq/scfastq/sample1_R2.fastq /home/john_doe/mbd1rnaseq/mm9.rsem sample1_R1.fastq sample1_R2.fastq
rsem-calculate-expression -p 8 --forward-prob 0.95 --bowtie2 --bowtie2-path /mnt/software/bowtie2/default --output-genome-bam --paired-end /home/john_doe/mbd1rnaseq/scfastq/sample2_R1.fastq /home/john_doe/mbd1rnaseq/scfastq/sample2_R2.fastq /home/john_doe/mbd1rnaseq/mm9.rsem sample2_R1.fastq sample2_R2.fastq
• 0 views • link
written by cpad0112 cpad0112

with gnu-parallel:

parallel --dry-run rsem-calculate-expression -p 8 --forward-prob 0.95 --bowtie2 --bowtie2-path /mnt/software/bowtie2/default --output-genome-bam   --paired-end ~/mbd1rnaseq/scfastq/{} ~/mbd1rnaseq/scfastq/{=s/R1/R2/=} ~/mbd1rnaseq/mm9.rsem {} {=s/R1/R2/=} ::: *R1.fastq


rsem-calculate-expression -p 8 --forward-prob 0.95 --bowtie2 --bowtie2-path /mnt/software/bowtie2/default --output-genome-bam --paired-end /home/john_doe/mbd1rnaseq/scfastq/sample1_R1.fastq /home/john_doe/mbd1rnaseq/scfastq/sample1_R2.fastq /home/john_doe/mbd1rnaseq/mm9.rsem sample1_R1.fastq sample1_R2.fastq
rsem-calculate-expression -p 8 --forward-prob 0.95 --bowtie2 --bowtie2-path /mnt/software/bowtie2/default --output-genome-bam --paired-end /home/john_doe/mbd1rnaseq/scfastq/sample2_R1.fastq /home/john_doe/mbd1rnaseq/scfastq/sample2_R2.fastq /home/john_doe/mbd1rnaseq/mm9.rsem sample2_R1.fastq sample2_R2.fastq
• 0 views • link
written by cpad0112 cpad0112

You'd need to change the initial glob, as your example needs the fastq files to be present in both $PWD and ~/mbd1rnaseq/scfastq/

• 0 views • link
written by Ram Ram

0 answers

No answers yet.

Log in to answer this question.

More posts like this

No similar posts found.