RNAseq for DE purpose
0
0
Entering edit mode
6 weeks ago
m.habib • 0

Hi all,

I am totally new in the bioinformatic analysis. I am working on a project that looks at DGE among different time treatments. Besides, there is no reference genome (meaning that I need a de novo assembly step).

So far, after struggling and navigating many website, I started using my terminal on MacBook and was able to do FastQC for the 118 reads (PE), I then followed this with MultiQC and both worked fine for me using the terminal. I am trying to go one step at a time. So I wanted to use trimmomatic to trim the reads but I failed!

I used the following command for one pair of reads and it worked but I do not know how to do this for the whole reads!

java -jar trimmomatic-0.39.jar PE -phred33 $i\1_R1_001.fastq.gz $i\1_R2_001.fastq.gz $i\1_R1_paired.fq.gz $i\1_R1_unpaired.fq.gz $i\1_R2_paired.fq.gz $i\1_R2_unpaired.fq.gz ILLUMINACLIP:contams_forward_rev.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 

I tried adding * but it did not work.

java -jar trimmomatic-0.39.jar PE -phred33 -threads 4 /Users/mohamedhabib/Desktop/RNAseq/*_R1_001.fastq.gz *_R1_001.fastq.gz out_*_R1_001_P.fastq.gz out_*_R1_001_U.fastq.gz out_*_R2_001_P.fastq.gz *_R2_001_U.fastq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:40:15

I will appreciate your help and advice on this specific topic and generally for further steps of analysis. Is there a template script (containing all the steps) that I can modify and use for my analysis?

Thanks

RNAseq • 812 views
ADD COMMENT
0
Entering edit mode

There is something weird about your commands: what are the \s doing after the variable names?

java -jar trimmomatic-0.39.jar PE -phred33 $i\1_R1_001.fastq.gz $i\1_...
#                                            ^                    ^
ADD REPLY
0
Entering edit mode

I just copied it as it is from another thread and modified it to my own data

ADD REPLY
0
Entering edit mode

What thread did you copy it from, what changes did you make and what was the logic behind each of those changes?

ADD REPLY
0
0
Entering edit mode

That answers 1 out of the 3 questions I asked.

ADD REPLY
0
Entering edit mode

I basically changed the file names to make sure this command is working. I picked one forward and reverse read pair and processed it as below: it worked but when I tried to trim all the samples (I thought maybe by adding a * to catch all the files!!) but it failed.

java -jar trimmomatic-0.39.jar PE -phred33 $i\1_R1_001.fastq.gz $i\1_R2_001.fastq.gz $i\1_R1_paired.fq.gz $i\1_R1_unpaired.fq.gz $i\1_R2_paired.fq.gz $i\1_R2_unpaired.fq.gz ILLUMINACLIP:contams_forward_rev.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

I hope that make sense because I am not an expert!

ADD REPLY
0
Entering edit mode

You need your command to run once per pair of files, not combine all pairs and then run just once. Read the thread you're trying to emulate and understand the loop part of it. Modify your code after you understand the loop.

ADD REPLY
0
Entering edit mode

I assume you are referring to 118 samlpes (and not reads). You can use a for loop to do this. See example here: Trimmomatic job script to run on multiple pair end read file

ADD REPLY
0
Entering edit mode

I used the following:

java -jar trimmomatic-0.39.jar PE -phred33 -threads 4 $i\_R1_001.fastq.gz $i\_R2_001.fastq.gz $i\_R1_001¬_paired.fastq.gz $i\_R1_001_unpaired.fastq.gz $i\_R2_001_paired.fastq.gz $i\_R2_001_unpaired.fastq.gz ILLUMINACLIP: TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 >> cmd_file

and I got:

Exception in thread "main" java.lang.ArrayIndexOutOfBoundsException: Index 1 out of bounds for length 1
ADD REPLY
1
Entering edit mode

Honestly, I think most RNASeq aligners are fairly robust to adapters. You could probably get away with not trimming at all

ADD REPLY
0
Entering edit mode

What names do your fastq files have? That first part of my answer will need to be adjusted to account for the differences in your file names so the for loop works right.

ADD REPLY
0
Entering edit mode

Here is an example: 19144R-03-01_S75_L004_R2_001.fastq.gz

ADD REPLY
0
Entering edit mode

Use following:

for i in `ls -1 *R1*.fastq.gz | sed 's/\_R1.fastq.gz//'` do echo trimmomatic PE -phred33 $i\_R1.fastq.gz $i\_R2.fastq.gz $i\_R1_paired.fq.gz $i\_R1_unpaired.fq.gz $i\_R2_paired.fq.gz $i\_R2_unpaired.fq.gz ILLUMINACLIP:contams_forward_rev.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 >> cmd_file; done

If the commands look correct then remove echo and >> cmd_file parts from the command to actually run the trimming.

ADD REPLY

Login before adding your answer.

Traffic: 2358 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6