SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
How to add a suffix to fastq file Wallysb01 Bioinformatics 9 08-27-2014 10:52 AM
how to remove 3'-adaptor sequence from illumina DGE expression data archory Bioinformatics 6 12-05-2011 07:55 AM
how to remove 3'-adaptor sequence from illumina DGE expression data archory Illumina/Solexa 0 11-29-2011 06:53 PM
Remove adapter sequence vini SOLiD 1 04-13-2011 10:28 AM
RNA-Seq: ABMapper: a suffix array-based tool for multi-location searching and splice- Newsbot! Literature Watch 0 12-21-2010 03:00 AM

Reply
 
Thread Tools
Old 03-13-2012, 04:10 AM   #1
alexd106
Junior Member
 
Location: Scotland

Join Date: Dec 2011
Posts: 4
Default remove suffix from fastq sequence ID

Dear all,

I have paired end illumina sequences in two large (20GiB) fastq files, one containing the forward reads, the other the reverse reads. Each file contains sequence IDs with either a /1 or /2 suffix. I would like to remove these suffixes (for some downstream analysis) from all reads and output 2 fastq files.

i.e.

change

@HWI-ST182_0249:5:1101:1093:2017#GTATGACG/1
NCAGCTGCAGGGAGTTAATTCACAGCAGTTGAGAGCCCTTGCTGTACCAACAAAGGGATGTGTGATCTCCCGGTCCCTCTGCCCCCTCCCCTCCCAGCCGC
+HWI-ST182_0249:5:1101:1093:2017#GTATGACG/1
BS\cacccegggehgghhhhh_ghhhhhhhhhhhhghhhhhhhhgghhhhhhhhhhhbghghghhhgeggedd`bb^bbbbbbaaaaaa_abaaabbaaaa

to

@HWI-ST182_0249:5:1101:1093:2017#GTATGACG
NCAGCTGCAGGGAGTTAATTCACAGCAGTTGAGAGCCCTTGCTGTACCAACAAAGGGATGTGTGATCTCCCGGTCCCTCTGCCCCCTCCCCTCCCAGCCGC
+HWI-ST182_0249:5:1101:1093:2017#GTATGACG
BS\cacccegggehgghhhhh_ghhhhhhhhhhhhghhhhhhhhgghhhhhhhhhhhbghghghhhgeggedd`bb^bbbbbbaaaaaa_abaaabbaaaa

I am new to bioinformatics and would appreciate a few pointers on the best way to get this done.
Thanks a million
Alex
alexd106 is offline   Reply With Quote
Old 03-13-2012, 05:12 AM   #2
rahularjun86
Member
 
Location: Frankfurt(M), Germany

Join Date: Jan 2011
Posts: 58
Default

Dear Alex,
You can use perl scripting, read the files, Split line if it is starting with @HWI or +HWI and print only the first part after splitting. And use else statement for printing rest of the sequence and quality lines as such.
Or you can use unix 'awk' set FS in the BEGIN and then print $1 part if line is starting with seq Id @HWI or +HWI.
Best wishes,
Rahul
__________________
Rahul Sharma,
Ph.D
Frankfurt am Main, Germany
rahularjun86 is offline   Reply With Quote
Old 03-13-2012, 06:56 AM   #3
ehlin
Member
 
Location: NYC

Join Date: Jan 2012
Posts: 12
Default

Quote:
Originally Posted by alexd106 View Post
Dear all,

I have paired end illumina sequences in two large (20GiB) fastq files, one containing the forward reads, the other the reverse reads. Each file contains sequence IDs with either a /1 or /2 suffix. I would like to remove these suffixes (for some downstream analysis) from all reads and output 2 fastq files.

i.e.

change

@HWI-ST182_0249:5:1101:1093:2017#GTATGACG/1
NCAGCTGCAGGGAGTTAATTCACAGCAGTTGAGAGCCCTTGCTGTACCAACAAAGGGATGTGTGATCTCCCGGTCCCTCTGCCCCCTCCCCTCCCAGCCGC
+HWI-ST182_0249:5:1101:1093:2017#GTATGACG/1
BS\cacccegggehgghhhhh_ghhhhhhhhhhhhghhhhhhhhgghhhhhhhhhhhbghghghhhgeggedd`bb^bbbbbbaaaaaa_abaaabbaaaa

to

@HWI-ST182_0249:5:1101:1093:2017#GTATGACG
NCAGCTGCAGGGAGTTAATTCACAGCAGTTGAGAGCCCTTGCTGTACCAACAAAGGGATGTGTGATCTCCCGGTCCCTCTGCCCCCTCCCCTCCCAGCCGC
+HWI-ST182_0249:5:1101:1093:2017#GTATGACG
BS\cacccegggehgghhhhh_ghhhhhhhhhhhhghhhhhhhhgghhhhhhhhhhhbghghghhhgeggedd`bb^bbbbbbaaaaaa_abaaabbaaaa

I am new to bioinformatics and would appreciate a few pointers on the best way to get this done.
Thanks a million
Alex
Hi Alex, while perl scripting is a good option, if you are new to bioinformatics there might be easier options for you. For example, FASTX-Toolkit:

http://hannonlab.cshl.edu/fastx_toolkit/
ehlin is offline   Reply With Quote
Old 03-13-2012, 07:00 AM   #4
alexd106
Junior Member
 
Location: Scotland

Join Date: Dec 2011
Posts: 4
Default

Hi Rahul,

Thank you very much for your suggestions. As i mentioned, I am new to bioinformatics and am just trying to teach myself some perl (and have never used awk). Would you mind providing a little more detail of the perl code you would use? No worries if not.

Cheers
Alex
alexd106 is offline   Reply With Quote
Old 03-13-2012, 07:18 AM   #5
kmcarr
Senior Member
 
Location: USA, Midwest

Join Date: May 2008
Posts: 1,172
Default

awk is good but sed might be faster and easier to learn.

Code:
sed -i.bak -e '/^[@+]HWI/ s/\/[12]$//' <yourFileName>
This sed script will look for lines starting with @HWI or +HWI, strip off either a /1 or /2 from the ends of those lines and save the result to the same file name as the original. The original file will be saved as <yourFileName>.bak.
kmcarr is offline   Reply With Quote
Old 03-13-2012, 07:52 AM   #6
alexd106
Junior Member
 
Location: Scotland

Join Date: Dec 2011
Posts: 4
Default

Thanks very much for the info.

All the best
Alex
alexd106 is offline   Reply With Quote
Old 03-13-2012, 08:01 AM   #7
rahularjun86
Member
 
Location: Frankfurt(M), Germany

Join Date: Jan 2011
Posts: 58
Default

Hi Alex,

Following is the perl code:
Code:
  1 use strict;
  2 use warnings;
  3 
  4 my $file_in=$ARGV[0];
  5 my $file_out=$ARGV[1];
  6 
  7 my $num=0;
  8 open I,"<$file_in" or die $!;
  9 open O,">$file_out" or die $!;
 10 
 11 do{
 12 
 13 my $f =<I>;
 14 chomp $f;
 15 
 16 if(($f =~ /^\@HWI/)||($f =~ /^\+HWI/))
 17      { $num++;
 18        my @s=split(/\//, $f);
 19        print O"$s[0]\n";
 20      }
 21 
 22 else
 23      {
 24        print O "$f\n";
 25         }
 26 
 27 }until eof(I);
 28 my $pr=$num/2;
 29 print "\nProcessed reads: $pr\n"
 30 
 31 
~                                                                                                                                                                    
~
Usage: perl program_name.pl Input_file.fq Out_file.fq
__________________
Rahul Sharma,
Ph.D
Frankfurt am Main, Germany

Last edited by rahularjun86; 03-13-2012 at 08:04 AM.
rahularjun86 is offline   Reply With Quote
Old 03-13-2012, 08:34 AM   #8
alexd106
Junior Member
 
Location: Scotland

Join Date: Dec 2011
Posts: 4
Default

Dear all, thanks for all the really useful suggestions. What a great community this is. I hope I can contribute sometime in the future when i have a little more experience.

[ehlin] I thought of using FASTX-Toolkit but couldn't see the appropriate tool. I looked at

$ fastx_renamer -h
usage: fastx_renamer [-n TYPE] [-h] [-z] [-v] [-i INFILE] [-o OUTFILE]
Part of FASTX Toolkit 0.0.10 by A. Gordon (gordon@cshl.edu)

[-n TYPE] = rename type:
SEQ - use the nucleotides sequence as the name.
COUNT - use simply counter as the name.

but it looks like the renaming is restricted to either a sequence or counter.

The sed and seemed to do the trick and I will look at the perl solution in an attempt the educate myself.
Cheers again
Alex
alexd106 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 08:35 PM.


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