SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
>90% aligned concordantly 0 times ChIP-seq Bowtie2 ge2sasag Illumina/Solexa 0 07-13-2017 04:04 AM
HISAT2 reads aligned concordantly exactly 1 time Beyounique Bioinformatics 2 05-18-2017 03:03 PM
aligned seq > 1 times and skipping reads because of seed mismatches carolW Bioinformatics 12 05-22-2014 11:48 PM
Bowtie2 - Concordantly tags. Coryza Bioinformatics 16 02-13-2014 06:00 AM

Reply
 
Thread Tools
Old 07-13-2017, 05:35 AM   #1
ge2sasag
Member
 
Location: Stockholm

Join Date: Apr 2013
Posts: 41
Question >90% aligned concordantly 0 times ChIP-seq Bowtie2

Hi,

I know this question have raised many times here and in other forums but I've tried everything other people suggested in previous posts and still can't figure it out where is the problem...... This is the output summary of Bowtie2 (I checked several times the target genome, hg19 and it's correct):

21404130 reads; of these:
21404130 (100.00%) were paired; of these:

21196512 (99.03%) aligned concordantly 0 times
104527 (0.49%) aligned concordantly exactly 1 time
103091 (0.48%) aligned concordantly >1 times
----
21196512 pairs aligned 0 times concordantly or discordantly; of these: 42393024 mates make up the pairs; of these:

906397 (2.14%) aligned 0 times
28296440 (66.75%) aligned exactly 1 time
13190187 (31.11%) aligned >1 times
97.88% overall alignment rate

I have paired-end ChIP-seq data, 50 bp reads. These are my steps (in Galaxy):

1) I groomed the fastq files to get fastqsanger (I checked if it's correct: Input FASTQ quality scores type --> Sanger & Illumina 1.8+)

2) FastQC is ok for all samples, some adapter contamination.

3) I trimmed using either TrimGalore or Trimmomatic (I get the same result using both) using forward and reverse files for each sample. I checked and I didn't missed up the forward and reverse samples (I run Bowtie2 using first file 1 and file 2 and also using first file 2 and second file 1 and I get the same result...)

4) I mapped the trimmed files to hg19 genome using Bowtie2 with option -fr and I get the result shown above. No matter what I try, I always get the same with the different samples that I have. All of them have around 32% of reads with adapters which I trimmed. I have no idea what to do with this.

Any suggestions or idea of where is the mistake? Am I missing something?

Thank you very much in advance.

Gema
ge2sasag is offline   Reply With Quote
Old 07-13-2017, 07:41 AM   #2
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,380
Default

It appears that your reads are not mapping concordantly. One possible reason is when you did the trimming you may not have used both files (R1/R2) at the same time which could have lead to reads being out of order in the two files. You should always trim paired-end data files together.

Last edited by GenoMax; 07-13-2017 at 07:44 AM.
GenoMax is offline   Reply With Quote
Old 07-13-2017, 07:47 AM   #3
ge2sasag
Member
 
Location: Stockholm

Join Date: Apr 2013
Posts: 41
Default

Hi GenoMax, thanks for your reply.

Both files were trimmed together with the option paired-end library... that's why I don't know the reason of disconcordant reads....
ge2sasag is offline   Reply With Quote
Old 07-13-2017, 07:48 AM   #4
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,380
Default

Have you examined the BAM files using a genome browser (e.g. IGV) to see if the two reads are mapping within expected distance (it appears that they perhaps are not, if the concordance message is to be believed).
GenoMax is offline   Reply With Quote
Old 07-13-2017, 07:55 AM   #5
ge2sasag
Member
 
Location: Stockholm

Join Date: Apr 2013
Posts: 41
Default

http://imgur.com/a/aFZ1h

this is an example of how the reads look like
ge2sasag is offline   Reply With Quote
Old 07-13-2017, 08:24 AM   #6
ge2sasag
Member
 
Location: Stockholm

Join Date: Apr 2013
Posts: 41
Default

Any idea??
ge2sasag is offline   Reply With Quote
Old 07-13-2017, 08:56 AM   #7
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,380
Default

It is difficult to know from that image since unlike IGV the two reads that belong to a fragment are not identifiable (IGV can display them as pairs). Default values for -I and -X are 0 and 500. Did you change either of them?

If you have the read files available locally can you check the insert sizes using BBMap program (it is pure Java and will run on PC/Mac)? See this answer to get the commands. You can add reads=1000000 to the command to use 1 million reads (instead of the full dataset).
GenoMax is offline   Reply With Quote
Old 07-13-2017, 09:21 AM   #8
ge2sasag
Member
 
Location: Stockholm

Join Date: Apr 2013
Posts: 41
Default

for -I, I used 0, for -X I tried 500 and 1000 with no difference.
ge2sasag is offline   Reply With Quote
Old 07-13-2017, 10:01 AM   #9
ge2sasag
Member
 
Location: Stockholm

Join Date: Apr 2013
Posts: 41
Default

This is an example of the first reads from TrimGalore output files 1 and 2 (as shown in Galaxy with the "eye" button):

File1:

@HWI-ST1018:141:H0HVBADXX:1:1101:1235:1969 1:N:0:CAGATC CCCANGATCTGTTCCACAGGAGATAAGCAGATCTTACTCCAGAGACCACTG + BB<f#0<b<fbffbfbbbff<b&lt;<ffffiibfffbf<fbbb<fb<ffiff7< p="">

@HWI-ST1018:141:H0HVBADXX:1:1101:1134:1970 1:N:0:CAGATC TGGCNTATGAAGTTCAGTGTTCTTTGGCTTGTTAGTCAGAACTGTTGC + BBBF#0BFFFFFFIFFFIIIIIFIIIFFIIIIIIIIIIIIIIFIIIII

@HWI-ST1018:141:H0HVBADXX:1:1101:1153:1984 1:N:0:CAGATC ACCTCTGTCTCCCAGGTTCAAGCGATTCTCCTGCCTCAGCCTCCCGAGT + BBBFFF<0B<b<fbfb&lt;0<ffb0fbb70bbbb0bff&lt;<bb<bfbfbfbb< p="">

File2:

@HWI-ST1018:141:H0HVBADXX:2:1101:1184:1971 1:N:0:CAGATC CTGATCAGAGGAGGAACATGACTAATCTATGGGCAGCCTACACTGAAGGC + BBBFFFFFFFFFFIIIIIIIIIIIIIFIIIIIIIIIIIIIIIIIIIIIII

@HWI-ST1018:141:H0HVBADXX:2:1101:1423:1966 1:N:0:CAGATC AATGGGTAGGTAAATGGATGGCTGGGTGAATGGATGGGTGGGTGGATTGGC + BBBFFFFFFFFFFIIIIFFIIFIIIIFFFIIIIFFFIIBFFIBFFFIIII7

@HWI-ST1018:141:H0HVBADXX:2:1101:1596:1990 1:N:0:CAGATC GATATCTTTTGTTTGTAGATATCTTTTCTAAGGCCCACATTCAGTGCAGAC + BBBFFBFFFFBBBBF<bfiiiiiiiiiiiifiifffiifiiiff0b<f7ff< p="">

Now I realize... it seems that after TrimGalore, the reads are not sorted?? TrimGalore is supposed to sort the reads right?
ge2sasag is offline   Reply With Quote
Old 07-13-2017, 10:12 AM   #10
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,380
Default

Something does not appear to be right. Did you get two output files (R1 trimmed and R2 trimmed) from the two (R1/R2 original) that you used as input for the trimming? Your file 2 does not contain any reads with 2 as read identifier (see example below).

e.g. @HWI-ST1018:141:H0HVBADXX:1:1101:1235:1969 1:N:0:CAGATC in File 1
should have a corresponding
@HWI-ST1018:141:H0HVBADXX:1:1101:1235:1969 2:N:0:CAGATC, in File 2. (note the bold read numbers).

Your example has no reads that show that property, if you are posting data from R1/R2 files for the same sample.

This is likely the reason why you are seeing discordant alignments.

You can either re-trim the data using a paired-end data aware aligner (you could use bbduk.sh from BBMap suite) or use "repair.sh" from BBMap to re-pair the trimmed files. Unfortunately you would need to do this outside galaxy.

Last edited by GenoMax; 07-13-2017 at 10:19 AM.
GenoMax is offline   Reply With Quote
Old 07-13-2017, 10:21 AM   #11
ge2sasag
Member
 
Location: Stockholm

Join Date: Apr 2013
Posts: 41
Default

But in file 2, there are no reads with 2:N:0 tag, the number 2 is here (in bold):

@HWI-ST1018:141:H0HVBADXX:2:1101:1184:1971 1:N:0:CAGATC CTGATCAGAGGAGGAACATGACTAATCTATGGGCAGCCTACACTGAAGGC + BBBFFFFFFFFFFIIIIIIIIIIIIIFIIIIIIIIIIIIIIIIIIIIIII

I already used a paired-end data aware aligner, Trim Galore, and those examples are the beginning of the two output files. Is it possible that the reads from the original FASTQ files are in the same orientation instead of -fr?
ge2sasag is offline   Reply With Quote
Old 07-13-2017, 10:30 AM   #12
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,380
Default

The 2 that you highlight in your example above refers to the lane number where this sample ran on a flowcell (see this).

Do you actually have paired-end data? i.e. Paired-end data should have two files per sample that have names like Sample_R1.fq.gz and Sample_R2.fq.gz. Can you post the names of your original files?

The two examples you posted above seem to contain data from one sample (or two independent samples) that ran in lanes 1 and 2. They are not paired-end data files.

Last edited by GenoMax; 07-13-2017 at 10:34 AM.
GenoMax is offline   Reply With Quote
Old 07-13-2017, 10:40 AM   #13
ge2sasag
Member
 
Location: Stockholm

Join Date: Apr 2013
Posts: 41
Default

aha! ok thanks for the clarification!

an example of original fastq file names:

INPUT_1_130506_BH0HVBADXX_P382_108_index8_1.fastq
INPUT_2_130506_BH0HVBADXX_P382_108_index8_1.fastq
ge2sasag is offline   Reply With Quote
Old 07-13-2017, 10:43 AM   #14
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,380
Default

I don't think these are paired-end data. These appear to be single-end data files for two samples (INPUT_1 and INPUT_2) and should be aligned independently of each other.
GenoMax is offline   Reply With Quote
Old 07-13-2017, 10:49 AM   #15
ge2sasag
Member
 
Location: Stockholm

Join Date: Apr 2013
Posts: 41
Default

Maybe that's the reason of the strange things.... I was confused by the numbers 1 and 2 in the sample name!! Thanks a lot for your help and time!!!
ge2sasag is offline   Reply With Quote
Old 07-13-2017, 10:57 AM   #16
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,380
Default

To be sure check with the experimental folks (if you did not do this experiment) to confirm sample provenance and how that correlates to the data files you have. I suggest that you go back and start the analysis process over. Trim and then align these files independently.
GenoMax is offline   Reply With Quote
Old 07-13-2017, 10:59 AM   #17
ge2sasag
Member
 
Location: Stockholm

Join Date: Apr 2013
Posts: 41
Default

sure, I will do it! I'm collaborating with another group and they just gave me the files...

I will re-do everything from the beginning and see how it goes.

THANKS!!
ge2sasag 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 06:51 PM.


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