Question: shell question: How to run multiple paired end alignment in shell?
0
gravatar for Jon
6 weeks ago by
Jon90
United States
Jon90 wrote:

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
ADD COMMENTlink modified 6 weeks ago • written 6 weeks ago by Jon90
1

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

ADD REPLYlink written 6 weeks ago by finswimmer11k
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
ADD REPLYlink modified 6 weeks ago • written 6 weeks ago by Jon90
1

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.

ADD REPLYlink modified 6 weeks ago • written 6 weeks ago by RamRS21k

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.
ADD REPLYlink modified 6 weeks ago by RamRS21k • written 6 weeks ago by Jon90

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.

ADD REPLYlink written 6 weeks ago by RamRS21k

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.
ADD REPLYlink written 6 weeks ago by Jon90
1

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.

ADD REPLYlink modified 6 weeks ago • written 6 weeks ago by RamRS21k
$ 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
ADD REPLYlink written 6 weeks ago by cpad011211k
1

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
ADD REPLYlink written 6 weeks ago by cpad011211k

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

ADD REPLYlink written 6 weeks ago by RamRS21k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1216 users visited in the last hour