SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Adapter trimming with cutadapt Kulvait Bioinformatics 3 05-05-2015 03:48 AM
cutadapt quality trimming cutoff algorithm bongbimit Bioinformatics 1 08-13-2014 08:36 AM
3' trimming using cutadapt changes property of 5'. Alex Lee RNA Sequencing 0 04-17-2014 10:28 PM
Adapters trimming: Cutadapt vs Trimmomatic MafaldaSF Bioinformatics 8 03-20-2014 06:16 AM
Using CASAVA versus cutadapt for adapter trimming id0 Bioinformatics 0 08-08-2013 09:39 AM

Reply
 
Thread Tools
Old 06-23-2015, 12:49 PM   #1
lianov
Junior Member
 
Location: USA

Join Date: Jun 2014
Posts: 5
Default cutadapt trimming for multiple files

Hi,

First, I'm sorry if I missed this but I couldn't find another thread about this issue:

I would like to run cutadapt on over 200 different .fastq files, so I clearly need to automate this process but I am not sure how to do this.

What is the best way to tell cutadapt to remove the adaptors on all files within a directory and create an output for each of them? The manual only mentions running two files at one (if paired end sequencing was done)

In short, I need to find all "*.fastq", trim sequences and output as trim.fastq

Thanks!
lianov is offline   Reply With Quote
Old 06-23-2015, 02:00 PM   #2
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,283
Default

You can use one of the shell script ideas here to iterate over the set of files you have: http://stackoverflow.com/questions/1...-files-in-unix
and http://stackoverflow.com/questions/1...s-in-directory

If you are submitting to a job scheduler on a cluster a similar loop can be used.
GenoMax is offline   Reply With Quote
Old 06-24-2015, 06:48 AM   #3
mgogol
Senior Member
 
Location: Kansas City

Join Date: Mar 2008
Posts: 193
Default

I bet you could do something snazzy with gnu parallel.

http://figshare.com/articles/GNU_par...otebook/822138
mgogol is offline   Reply With Quote
Old 08-20-2015, 04:13 AM   #4
slefevre
Junior Member
 
Location: Norway

Join Date: Jun 2015
Posts: 3
Default

Hi - for starters, I am a completely newbie at this, so bare with me if it's basic.

I understand how to do a for loop with one object, but if you have paired-end reads, trim-galore takes two files as input, and I just can't figure out how to make a for loop including a list of pairs.... Would be easy if trim-galore accepted wildcards in the two filenames, but I am not sure of it does? Or is it possible to specify two objects to use as input? i.e.
for read1 in /*/xxxxR1.fastq
for read2 in /*/xxxxR2.fastq
do trim-galore etc etc "$read1$ "$read2"
done

Any thoughts?
slefevre is offline   Reply With Quote
Old 08-20-2015, 05:39 AM   #5
Michael.Ante
Senior Member
 
Location: Vienna

Join Date: Oct 2011
Posts: 120
Default

You can deduce the read2 variable from read1:

Code:
for read1 in */xxxxR1.fastq; do read2=$(echo $read1| sed 's/R1.fastq/R2.fastq/'); trim-galore etc etc $read1 $read2 ; done
Michael.Ante is offline   Reply With Quote
Old 08-20-2015, 05:40 AM   #6
mgogol
Senior Member
 
Location: Kansas City

Join Date: Mar 2008
Posts: 193
Default

Ah. The way I usually handle stuff like this is just I have a text file listing all the files - in this case, in two columns. Then I have a perl script with a for loop going through each line in the text file and generating the commands I want to run. It's useful later on in R or whatever as well.
mgogol is offline   Reply With Quote
Old 08-20-2015, 05:46 AM   #7
slefevre
Junior Member
 
Location: Norway

Join Date: Jun 2015
Posts: 3
Default

Quote:
Originally Posted by Michael.Ante View Post
You can deduce the read2 variable from read1:

Code:
for read1 in */xxxxR1.fastq; do read2=$(echo $read1| sed 's/R1.fastq/R2.fastq/'); trim-galore etc etc $read1 $read2 ; done
That makes sense - and even better, it seems to be working. Thank you so much!
slefevre is offline   Reply With Quote
Old 08-20-2015, 07:44 AM   #8
Michael.Ante
Senior Member
 
Location: Vienna

Join Date: Oct 2011
Posts: 120
Default

Quote:
Originally Posted by mgogol View Post
Ah. The way I usually handle stuff like this is just I have a text file listing all the files - in this case, in two columns. Then I have a perl script with a for loop going through each line in the text file and generating the commands I want to run. It's useful later on in R or whatever as well.
In case you want to speed up a bit, you can use such a file as input for GNU parallel:

Code:
parallel -j 8 --colsep '\t'  trim-galore etc etc {1} {2} :::: read-pairs-location.tsv
Just adjust the number of threads (-j) according to your CPUs.
Michael.Ante is offline   Reply With Quote
Old 11-13-2015, 12:44 PM   #9
yaseen.ladak
Junior Member
 
Location: Cambridge

Join Date: Nov 2013
Posts: 7
Default

Hello Michael,

I am trying to do exactly as you suggested just need a bit of help with that code.

what is --colsep and when you say '\t' means that you read-pairs-location.tsv is tab separated if it was a CSV you would write --colsep ',' please correct me if I am wrong. Also it means that it will use 8 threads for each of the trim_galore command?

I am currently testing on a Quad code i7 on my mac that has hyper threading so in total I have 8 cores 4 physical and 4 virtual so I can use -j8 although I need to test for what value of thread I get the best performance as using maximum threads generally gives lower performance.

Thanks,
Yaseen
yaseen.ladak is offline   Reply With Quote
Old 11-13-2015, 01:02 PM   #10
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,283
Default

That --colsep option should work. Give it a try Otherwise convert your file into tsv format and then use the original \t option.

Having 8 cores does not mean that you will be able to use them efficiently. They are connected to your storage subsystem through a common bus which can only allow a certain amount of data to be read/written.

Again experiment with a small subset of data (test file) starting with 4 cores and go up and down in number to see what number works best for your setup.
GenoMax is offline   Reply With Quote
Old 11-13-2015, 02:45 PM   #11
yaseen.ladak
Junior Member
 
Location: Cambridge

Join Date: Nov 2013
Posts: 7
Default

Thank you GenoMax.

I ran the following command:

parallel -j 3 --colsep '\t' trim_galore --output_dir tg/3 --paired --fastqc -a CTGTCTCTTATA -a2 CTGTCTCTTATA {1} {2} :::: tg.tsv

Is there a way that I can give a name to this command? and how would i check the number of threads being used. I am using a mac and when I check the activity monitor i cannot see 3 i.e. number of threads. Had i named this command it would have shown up in the activity monitor and or using top? Indeed this command is running as it is generating output. The tg.tsv has 4 files names with their complete paths i.e. 2 paired end samples.

Can someone please advise? thanks
yaseen.ladak is offline   Reply With Quote
Old 11-13-2015, 02:52 PM   #12
yaseen.ladak
Junior Member
 
Location: Cambridge

Join Date: Nov 2013
Posts: 7
Default

In my activity monitor when i run the above command its shows gzip command with 1 thread as i think it uncompressed the fastq.gz file initially? Whats the best to see the number of threads run by this command on a mac or a linux. I have a mac.
Can I name this command to something so i can quickly find it in top etc?
yaseen.ladak is offline   Reply With Quote
Old 11-13-2015, 04:54 PM   #13
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,283
Default

When you run "top" processes that are actively consuming CPU cycles should show up at the top of the list. You can also use the activity monitor: https://support.apple.com/en-us/HT201464

Your file contains the two file name pairs separated by a tab on each line, right?
GenoMax is offline   Reply With Quote
Old 11-14-2015, 02:11 AM   #14
yaseen.ladak
Junior Member
 
Location: Cambridge

Join Date: Nov 2013
Posts: 7
Default

When I do a top or a activity monitor it shows gzip most of the time and the number of threads as 1 even when i run with 3 threads.

Yes tsv in a text editor, with two reads each separated by a tab. I did not put two samples as my mac gives me out of space error so for testing I can only use 1 sample with paired end read. I then labelled the extension of the file with .tsv The below is the section how my tab separated file looks like, the gap between the two files is a tab.

/Users/Yaseen/S1R1.fastq.gz /Users/Yaseen/S1R2.fastq.gz
yaseen.ladak is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 03:27 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2017, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO