Question: Remove duplicates reads Ids
2
gravatar for BioGeek
2.1 years ago by
BioGeek140
BioGeek140 wrote:

I mapped reads with

bwa mem -M -t 40 allCombinedFinalSet.fa Seq.R1.fastq Seq.R2.fastq > aln.sam

Extracted the mapped reads

samtools view -f 0x2 -b aln.bam > output.bam

Extracted the fastq

bamToFastq -i output.bam -fq R1.fq -fq2 R2.fq 

grep @HISEQ578:1035:HJ2KCBCXX:1:1104:14672:39678/1 R1.fq             []
@HISEQ578:1035:HJ2KCBCXX:1:1104:14672:39678/1
@HISEQ578:1035:HJ2KCBCXX:1:1104:14672:39678/1
@HISEQ578:1035:HJ2KCBCXX:1:1104:14672:39678/1

I notice it has duplicated ....

I think this because read was mapped twice (i.e. BWAmem).

I tried fastuniq but it does not remove the duplicated reads.

Can you please help me to remove duplicated reads from fastq files.

ids mapping fastq duplicates reads • 1.5k views
ADD COMMENTlink modified 2.1 years ago • written 2.1 years ago by BioGeek140

What exactly are you trying to do?

ADD REPLYlink written 2.1 years ago by genomax67k

I am trying to mapped the allReads against "filtered" genome and extract mapped fastq for re-assembly.

ADD REPLYlink written 2.1 years ago by BioGeek140
0
gravatar for Santosh Anand
2.1 years ago by
Santosh Anand4.7k
Santosh Anand4.7k wrote:

Those duplicates are secondary/supplementary alignments. You may try Picard SamToFastq for a better control over what to convert in Fastq from BAM. See the discussion here: http://bedtools.readthedocs.io/en/latest/content/tools/bamtofastq.html

ADD COMMENTlink written 2.1 years ago by Santosh Anand4.7k
0
gravatar for k.balaji93
2.1 years ago by
k.balaji9310
Penn State
k.balaji9310 wrote:

Read Samtools man for flag descriptions. You want to filter out secondary and chimeric alignments, using samtools in the bam file. Check out Picard explain flag to get the exact flag to keep only primary alignments.

ADD COMMENTlink modified 2.1 years ago • written 2.1 years ago by k.balaji9310
0
gravatar for Pierre Lindenbaum
2.1 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum120k wrote:

(assuming that the reads in R1.fq are in the same order + same number than R2.fq)

paste <(cat R1.fq | paste - - - - )  <(cat R2.fq | paste - - - - ) | LC_ALL=C sort -t $'\t' -k1,1  | tr "\t" "\n" > interleaved.fq
ADD COMMENTlink modified 2.1 years ago • written 2.1 years ago by Pierre Lindenbaum120k
bio@bio214b[biogeek] paste <(cat see.R1.fq | paste - - - - )  <(cat see.R2.fq | paste - - - - ) | LC_ALL=C sort -t $'\t' -k1,1  | tr "\t" "\n" > interleaved.fq
bio@bio214b[biogeek] grep HJ2KCBCXX:1:1104:14672:39678/1 interleaved.fq
@HISEQ578:1035:HJ2KCBCXX:1:1104:14672:39678/1
@HISEQ578:1035:HJ2KCBCXX:1:1104:14672:39678/1
@HISEQ578:1035:HJ2KCBCXX:1:1104:14672:39678/1
bio@bio214b[biogeek] grep HJ2KCBCXX:1:1104:14672:39678/2 interleaved.fq
@HISEQ578:1035:HJ2KCBCXX:1:1104:14672:39678/2
@HISEQ578:1035:HJ2KCBCXX:1:1104:14672:39678/2
@HISEQ578:1035:HJ2KCBCXX:1:1104:14672:39678/2

The duplicated still exists.

ADD REPLYlink modified 2.1 years ago • written 2.1 years ago by BioGeek140

If all you need to do is remove duplicates then use dedupe.sh from BBMap suite. dedupe.sh in=orig.fq out=dedupe.fq Should work with in1= in2= out1= out2= for PE files.

ADD REPLYlink modified 2.1 years ago • written 2.1 years ago by genomax67k
$ dedupe.sh in1=see.R1.fq in2=see.R2.fq out1=tt.R1.fq out2=tt.R2.fq          []
java -Djava.library.path=/home/Tools/bbmap/jni/ -ea -Xmx200291m -Xms200291m -cp /home/Tools/bbmap/current/ jgi.Dedupe in1=see.R1.fq in2=see.R2.fq out1=tt.R1.fq out2=tt.R2.fq
Executing jgi.Dedupe [in1=see.R1.fq, in2=see.R2.fq, out1=tt.R1.fq, out2=tt.R2.fq]

Exception in thread "main" java.lang.RuntimeException: Unknown parameter out1=tt.R1.fq
    at jgi.Dedupe.<init>(Dedupe.java:392)
    at jgi.Dedupe.main(Dedupe.java:82)

It gave an error !!!

ADD REPLYlink modified 2.1 years ago • written 2.1 years ago by BioGeek140

While dedupe.sh will work to remove duplicate read ID's I realized that it will have an additional unintended consequence. It will actually de-duplicate your input dataset based on sequence. If you are able to accept that caveat then read on.

dedupe.sh requires input to be interleaved in case of paired-end reads so you would be able to use it like this:

  1. Interleave reads
  2. de-dupe reads
  3. Split the deduped reads into separate R1/R2 files

reformat.sh in1=R1.fq in2=R2.fq out=stdout.fq | dedupe.sh in=stdin.fq out=stdout.fq | reformat.sh in=stdin.fq out1=newR1.fq out2=newR2.fq

I am going to tag Brian Bushnell to see if he has any other suggestion.

My apologies for not thinking this through completely.

ADD REPLYlink modified 2.1 years ago • written 2.1 years ago by genomax67k

I don't really have a tool for this purpose, once you already have a bam file. Instead, I'd recommend using BBMap, like this:

bbmap.sh t=40 ref=allCombinedFinalSet.fa in1=Seq.R1.fastq in2=Seq.R2.fastq outm=mapped.fq outu=unmapped.fq

The mapped reads will go to mapped.fq, and the unmapped reads will go to unmapped.fq. Both files will be interleaved, and pairs will be kept together, so if only one read in a pair maps to the reference, both reads will still go to mapped.fq.

ADD REPLYlink written 2.1 years ago by Brian Bushnell16k
0
gravatar for BioGeek
2.1 years ago by
BioGeek140
BioGeek140 wrote:
    #!/usr/bin/perl-w
use strict;
use warnings;

my %id2seq=();
my $key = '';
$/ = '@HISEQ';
my $add='@HISEQ';
while(<>){
    chomp;
    my ($defLine, @seqLines) = split /\n/, $_;
    my $sequence = join("\n",@seqLines);
        $id2seq{$defLine} = $sequence;
}

foreach(keys %id2seq){
    my $name="$add"."$_";  
    print join("\n",$name,$id2seq{$_}),"\n";
}

I tried my above mentioned perl script to remove, and in first look it seems to be working.

bio@bio214b[biogeek] grep HJ2KCBCXX:1:1104:14672:39678/2 test.R2.fq 
@HISEQ578:1035:HJ2KCBCXX:1:1104:14672:39678/2
bio@bio214b[biogeek] grep HJ2KCBCXX:1:1104:14672:39678/1 test.R1.fq
@HISEQ578:1035:HJ2KCBCXX:1:1104:14672:39678/1

Is that safe to use? It seems distorting the format.

ADD COMMENTlink modified 2.1 years ago • written 2.1 years ago by BioGeek140

It seems distorting the format

What exactly is it doing?

ADD REPLYlink written 2.1 years ago by genomax67k
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: 996 users visited in the last hour