1
$\begingroup$

I have numerous gzip-compressed paired-end Fastq files with ChIP-seq data that I would like to align with bowtie2. I confirmed that bowtie2 takes .gz files, although I do see this warning when I run for two samples:

Warning: gzbuffer added in zlib v1.2.3.5. Unable to change buffer size from default of 8192.

My files are named like so.

10_S10_L002_R1_001.fastq.gz  12_S12_L002_R1_001.fastq.gz  2_S2_L001_R1_001.fastq.gz  4_S4_L001_R1_001.fastq-003.gz  6_S6_L001_R1_001.fastq-002.gz  8_S8_L002_R1_001.fastq.gz  
10_S10_L002_R2_001.fastq.gz  12_S12_L002_R2_001.fastq.gz  2_S2_L001_R2_001.fastq.gz  4_S4_L001_R2_001.fastq-002.gz  6_S6_L001_R2_001.fastq-003.gz  8_S8_L002_R2_001.fastq.gz  
11_S11_L002_R1_001.fastq.gz  1_S1_L001_R1_001.fastq.gz    3_S3_L001_R1_001.fastq.gz  5_S5_L001_R1_001.fastq-002.gz  7_S7_L002_R1_001.fastq.gz      9_S9_L002_R1_001.fastq.gz
11_S11_L002_R2_001.fastq.gz  1_S1_L001_R2_001.fastq.gz    3_S3_L001_R2_001.fastq.gz  5_S5_L001_R2_001.fastq-003.gz  7_S7_L002_R2_001.fastq.gz      9_S9_L002_R2_001.fastq.gz

The bash commands I tried are below. When I run this test command:

for i in $(ls *.fastq*.gz | grep -v 6 | rev | cut -c 16- | rev | uniq); do 
    echo $i; 
done

...I get this output.

12_S12_L002_R1_001.fastq.gz  12_S12_L002_R2_001.fastq.gz 

When I run this:

for i in $(ls *.fastq*.gz | rev | cut -c 16- | rev | uniq); do 
    echo $i 
    bowtie2 -q -p20 -x /run/media/punit/data1/BowtieIndex/hg38 -1 ${i}R1_001.fastq.gz -2 ${i}R2_001.fastq.gz | samtools view -bS - >${i%}.bam
done

I do see one BAM for my paired end files. But as I see the second set of files which are named:

6_S6_L001_R1_001.fastq-002.gz  6_S6_L001_R2_001.fastq-003.gz

if I run the same script over them, I would get an error as the loop can't read the paired end file.

So the tedious way is I can make a separate folder for these sets of data with fastq-002.gz and fastq-003.gz and can run the script by tweaking the cut parameter.

But I would like to know what is the programmatic way of resolving the issue.

$\endgroup$
6
  • 1
    $\begingroup$ I strongly encourage you to use something like snakeMake rather than making a shell script. It's both much easier to ensure the proper files are paired and it scales to use in clusters and the cloud. $\endgroup$
    – Devon Ryan
    Mar 25 '19 at 14:24
  • $\begingroup$ Can't you just rename the files? $\endgroup$
    – terdon
    Mar 25 '19 at 15:06
  • $\begingroup$ @terdon that's how our collaborator gave us and i dont want to mess up with the file name yes i can if i dont find any other way.. $\endgroup$
    – kcm
    Mar 25 '19 at 15:13
  • $\begingroup$ " strongly encourage you to use something like snakeMake rather than making a shell script. It's both much easier to ensure the proper files are paired and it scales to use in clusters and the cloud." i read a bit about it never tried it , i need to learn and give it a try $\endgroup$
    – kcm
    Mar 25 '19 at 15:14
  • 1
    $\begingroup$ It's good practice to treat files from collaborators as read-only, so not changing the filenames is a good idea. But if the files are inconsistently named, it's a very good idea to create consistently named symlinks to those files. Well organized filenames and directory structures will make it SO NICE AND EASY to automate batch data processing. $\endgroup$ Mar 26 '19 at 18:31
3
$\begingroup$

A very simple approach, but one which will work with the names you show, is to remove any occurrence of 00N (where N is any number) from right before the .gz. To do so, edit your script to:

for i in $(ls *.fastq*.gz | sed 's/00[0-9]\.gz/.gz/' | rev | cut -c 16- | rev | uniq); do 
  bowtie2 -q -p20 -x /run/media/punit/data1/BowtieIndex/hg38 \
    -1 ${i}R1_001.fastq.gz \
    -2 ${i}R2_001.fastq.gz | 
        samtools view -bS - >${i%}.bam"
done
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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