SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
FastX-toolkit liu_xt005 Bioinformatics 13 10-11-2014 04:52 AM
FASTX-Toolkit: quality score value thinkRNA Bioinformatics 13 09-30-2014 09:25 AM
Has anyone used FastX Toolkit with IonTorrent data? gkijak Ion Torrent 7 06-05-2013 05:16 PM
Filtering Illumina GAIIx reads using fastx rubi Bioinformatics 8 05-02-2011 03:54 AM
matching unmapped paired SOLiD reads smarkel General 6 10-09-2009 11:49 AM

Reply
 
Thread Tools
Old 04-24-2014, 11:51 AM   #41
ericaramos
Junior Member
 
Location: SP, Brazil

Join Date: Apr 2014
Posts: 4
Default

Quote:
Originally Posted by carmeyeii View Post
Dear btmb,
I'm afraid I still cannot run it. Sorry to keep bothering?

I have corrected tabs and spaces to avoid getting the Unexpected indent Error,

but now I get:



Thanks again for any help,

Carmen
Quote:
Originally Posted by SES View Post
If you look through the discussion above you can see that a number of people had similar issues, and this script doesn't appear to be maintained. I think the best solution may be to find another approach unless you want to work on that shell/python code.

Did you try the tool Pairfq that was mentioned in the thread above? I'd be happy to help with this if you run into any issues. We can help with the other approach as well, but it is hard to see what the issue is and it's also a challenge to keep code updated on a forum such as this.



...................................................................................................................................Ok, I didn't try using Pairfq, but I will.

Thank you for the answer!
ericaramos is offline   Reply With Quote
Old 04-25-2014, 11:14 AM   #42
ericaramos
Junior Member
 
Location: SP, Brazil

Join Date: Apr 2014
Posts: 4
Default

Quote:
Originally Posted by SES View Post
If you look through the discussion above you can see that a number of people had similar issues, and this script doesn't appear to be maintained. I think the best solution may be to find another approach unless you want to work on that shell/python code.

Did you try the tool Pairfq that was mentioned in the thread above? I'd be happy to help with this if you run into any issues. We can help with the other approach as well, but it is hard to see what the issue is and it's also a challenge to keep code updated on a forum such as this.
--------------------------------------------------------------------------------------------------------------------------------------
Pairfq worked pretty well!! Thank you!
ericaramos is offline   Reply With Quote
Old 05-28-2014, 02:46 AM   #43
ranu1
Junior Member
 
Location: Delhi

Join Date: May 2014
Posts: 2
Default

After removing the adapters from cutadapt i got unsymmetrical pair end file so I want to know the script that could remove the orphan reads and make the data symmetric although I made it using hash but its very slow.The above mention script is showing error..
ranu1 is offline   Reply With Quote
Old 05-28-2014, 05:37 AM   #44
SES
Senior Member
 
Location: Vancouver, BC

Join Date: Mar 2010
Posts: 275
Default

Quote:
Originally Posted by ranu1 View Post
After removing the adapters from cutadapt i got unsymmetrical pair end file so I want to know the script that could remove the orphan reads and make the data symmetric although I made it using hash but its very slow.The above mention script is showing error..
We will need some more details in order to help. For example, which script are you referring to, the Python script mentioned on the first page of this thread? If that is the script you are attempting to use, I don't think you'll be able to get it working without some code changes, as mentioned above.

Also, what do you mean when you say the script is showing error? It is not possible to know what the issue is based on that information alone.
SES is offline   Reply With Quote
Old 05-28-2014, 07:44 AM   #45
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

BBTools has a tool to quickly re-pair arbitrarily disordered reads based on their names.

For interleaved reads:

repair.sh in=reads.fq out=fixed.fq outsingle=single.fq

For paired reads in two files:

repair.sh in1=read1.fq in2=read2.fq out1=fixed1.fq out2=fixed2.fq outsingle=single.fq

You can also repair simple broken interleaving much faster and with less memory, but this will not fix arbitrarily disordered reads, just reads that were interleaved and had some of the reads thrown away:

bbsplitpairs.sh in=reads.fq out=fixed.fq outsingle=single.fq fixinterleaving

Last edited by Brian Bushnell; 02-13-2015 at 09:31 AM.
Brian Bushnell is offline   Reply With Quote
Old 05-30-2014, 02:50 AM   #46
SES
Senior Member
 
Location: Vancouver, BC

Join Date: Mar 2010
Posts: 275
Default

Quote:
Originally Posted by Brian Bushnell View Post
BBTools now has a tool to quickly re-pair arbitrarily disordered reads based on their names.

For interleaved reads:

repair.sh -Xmx29g in=reads.fq out=fixed.fq outsingle=single.fq

For paired reads in two files:

repair.sh -Xmx29g in1=read1.fq in2=read2.fq out1=fixed1.fq out2=fixed2.fq outsingle=single.fq

You can also repair simple broken interleaving much faster and with less memory, but this will not fix arbitrarily disordered reads, just reads that were interleaved and had some of the reads thrown away:

bbsplitpairs.sh in=reads.fq out=fixed.fq outsingle=single.fq fixinterleaving
I tried to test this latest script but I can't get it to run. There seems to be a bug. I get: "./repair.sh: line 73: splitpairs: command not found." Beyond that, it seems like it has to be run in place and can't be moved because it runs other shell scripts. I noticed that the script runs ulimit every time it is executed, even if you just want a menu. This is kind of annoying on a cluster because it calculates the number of files for every mounted file system, which could take a while.

For reference, I was simply trying to compare this to Pairfq (a tool for pairing reads, which was mentioned above), since I am the author of that tool.
SES is offline   Reply With Quote
Old 05-30-2014, 09:01 AM   #47
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

SES,

Thanks for reporting that! I'll have a fixed version of the shellscript up later today. If you want to test it in the meantime, you can replace
repair.sh
with
java -ea -Xmx8g -cp /path/to/bbmap/current/ jgi.SplitPairsAndSingles

...where "-Xmx8g" can be adjusted to be around 80% of your node's physical memory.

For a menu, "repair.sh" or "repair.sh -h" will still work and they should not run ulimit.

But I'm a bit confused about your comment on calculating files in mounted file systems. My cluster is attached to several multi-petabyte file systems containing hundreds of millions of files, a couple of which are very slow, and "ulimit -v" always returns an answer instantly even if "ls" takes several seconds to respond.
Brian Bushnell is offline   Reply With Quote
Old 05-30-2014, 09:10 AM   #48
SES
Senior Member
 
Location: Vancouver, BC

Join Date: Mar 2010
Posts: 275
Default

Quote:
Originally Posted by Brian Bushnell View Post
SES,

Thanks for reporting that! I'll have a fixed version of the shellscript up later today. If you want to test it in the meantime, you can replace
repair.sh
with
java -ea -Xmx8g -cp /path/to/bbmap/current/ jgi.SplitPairsAndSingles

...where "-Xmx8g" can be adjusted to be around 80% of your node's physical memory.

For a menu, "repair.sh" or "repair.sh -h" will still work and they should not run ulimit.

But I'm a bit confused about your comment on calculating files in mounted file systems. My cluster is attached to several multi-petabyte file systems containing hundreds of millions of files, a couple of which are very slow, and "ulimit -v" always returns an answer instantly even if "ls" takes several seconds to respond.
Thanks for the response, I'll give this a try when I have a chance. About ulimit, you are right, that command does return instantly. It must be another command that is printing the number of files for each file system when I run the script but I can't figure out what that is right now.
SES is offline   Reply With Quote
Old 06-03-2014, 12:31 PM   #49
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

SES,

I uploaded v32.27 which fixes the bug you found in repair.sh. I'd be interested in knowing whether it still prints directory information on your system.

-Brian
Brian Bushnell is offline   Reply With Quote
Old 01-04-2015, 08:49 PM   #50
luqmanaslam
Junior Member
 
Location: USA

Join Date: Oct 2014
Posts: 2
Default Pairfq

Hi! I am having issue in using pairfq. Probably, it does not identify Illumina 1.9 fastq format and gives me following error.
ERROR: Could not determine FASTA/Q format. Please see https://github.com/sestaton/Pairfq or the README for supported formats. Exiting.

could you please update it to solve this problem?
If someone has and can share tool and/or script to sort and separate paired-end common and unique reads.

My reads and header identifiers look like this
Read1.fastq
@6_1101_1236_2393_1
AATTCATAAAAAACAACTGAATTGGATCACTGTAGTACAAATTGCAAACATTCACTATGGGCAAACAAGNNNNNNNNNNNNNNNNNNNNNNNNNN
+
FDDHHHHFJJJIIJJIIIJGHHIJJEHHJJIJHIGIJJJJJHIIIJJJIJIJGIIJJJJIJHIIJJGHH##########################
@6_1101_1728_2439_1
AATTCATCGGGAACGACCTGTTGCCTTTTCAAATCTTTCTAAGTTATTCTGTTCAGCTGCAAAGACTTGACATATTTNNNNNNNNNNNNNNNNNN
+
FFFHHGHHJJJGIJJJJJJJJJJJJJJJFGGIJJIIIIJJJIJIIGJJJJJJJJJJIIJJJFJGFHDHHFFFFCEEE##################

Read2.fastq
@6_2316_20850_100973_2
AATAAAATCAATCCAGACTAGCGGCACATTTTGACTGTTAAGTTTGAACTTCCTAAAATCTGTGCAGGCTTTTAAGCTGTACTGTTTTTCTT
+
HFHHHIJJIJJIJJIDHIJJFHIBHGJIBHIJICHIIJJ>CBFBFH=CFIIIJI).@DHHGEEEHDFBFECEEC@>@CDCCDBCBCDDDDDD
@6_2316_21147_100977_2
GCTGTTTTTAAATTTAAAATTTAAAAAGTGCTTTTTGTTGTTGTTTTTTTAAAAGGAAAAGGCTTCCTATAACTTCTCATGCTGGACAACAC
+
HHHHHJJJJEGHJIJJJJJJJJIIHHIJHGHIJJJJJIJJJJGHIIJJJJJJEH=AHFF>DFDEDEEDDD>CDDEDDDDAACDCCBDAAABD


Thanks
luqmanaslam is offline   Reply With Quote
Old 01-04-2015, 10:41 PM   #51
SES
Senior Member
 
Location: Vancouver, BC

Join Date: Mar 2010
Posts: 275
Default

The expected pair information is missing from these reads, at least I have not seen data of this format. The solution is to add the pair information and then pair the reads.

Code:
pairfq addinfo -i Read1.fastq -o Read1_pairinfo.fastq -p 1
pairfq addinfo -i Read2.fastq -o Read2_pairinfo.fastq -p 2
Now you can use 'pairfq makepairs' to pair up the reads.

Code:
pairfq makepairs -f Read1_pairinfo.fastq \
-r Read2_pairinfo.fastq \
-fp Read1_p.fastq \
-rp Read2_p.fastq \
-fs Read1_s.fastq \
-rs Read2_s.fastq \
--stats
Hope that helps. Also, I'd be interested to know if these read IDs were generated from an Illumina instrument or a custom pipeline. If it is the former then I can add support for this type of data.

Last edited by SES; 01-05-2015 at 10:13 AM. Reason: fix typo
SES is offline   Reply With Quote
Old 01-05-2015, 06:17 PM   #52
luqmanaslam
Junior Member
 
Location: USA

Join Date: Oct 2014
Posts: 2
Default

Thanks for this, yes these IDs are not generated by any sequencing platform. These come out of Stacks pipeline when you use demultiplexing and quality filtering of RAD data. It will be good to consider this in the script which will widen the usage of program straight away for many other people like me.

Actually, I got another different issue. I started using a Linux machine from a different lab and I tried to submit script using qsub command and it says
-bash: qsub: command not found
when I write
which qsub
I get following
/usr/bin/which: no qsub in (/usr/lib64/qt-3.3/bin:/usr/local/bin:/bin:/usr/bin:/usr/local/sbin:/usr/sbin:/sbin:/home/user/bin)

When I tried
man qsub
it opens qsub help and instruction.
I don't understand why I am not able to submit my shell script
I am trying to use following codes for submission
qsub -q all.q -cwd test.sh

I will appreciate any help or guidance in this regard.
Thanks,
luqmanaslam is offline   Reply With Quote
Old 01-05-2015, 06:38 PM   #53
SES
Senior Member
 
Location: Vancouver, BC

Join Date: Mar 2010
Posts: 275
Default

Quote:
Originally Posted by luqmanaslam View Post
Thanks for this, yes these IDs are not generated by any sequencing platform. These come out of Stacks pipeline when you use demultiplexing and quality filtering of RAD data.
That is good to know, it would be helpful for them to follow the FASTQ specifications though.


Quote:
Actually, I got another different issue. I started using a Linux machine from a different lab and I tried to submit script using qsub command and it says
-bash: qsub: command not found
when I write
which qsub
I get following
/usr/bin/which: no qsub in (/usr/lib64/qt-3.3/bin:/usr/local/bin:/bin:/usr/bin:/usr/local/sbin:/usr/sbin:/sbin:/home/user/bin)

When I tried
man qsub
it opens qsub help and instruction.
I don't understand why I am not able to submit my shell script
I am trying to use following codes for submission
qsub -q all.q -cwd test.sh

I will appreciate any help or guidance in this regard.
Thanks,
Are you sure the cluster is using SGE as a queuing system, that you have permissions on the queue called "all.q" and that this queue exists? Hopefully, that will help you start to figure it out, but the best solution would be to talk to someone informed about the cluster, that will probably save you and others a lot of time. Good luck!
SES is offline   Reply With Quote
Old 01-05-2015, 06:39 PM   #54
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

qsub is just for clusters using UGE/SGE (as far as I know). If you are trying to operate on a cluster with a different scheduler, you would need a different command; but if you just want to run in Linux without submitting a job to a scheduler, then your command in this case would simply be "sh test.sh".
Brian Bushnell is offline   Reply With Quote
Old 07-30-2015, 03:13 AM   #55
Armin_P
Junior Member
 
Location: Austria

Join Date: Sep 2014
Posts: 4
Default Perl script to handle to issue

Since this has generally remained a common problem, I wrote a Perl script that removes the unpaired reads and matches the paired ones. Just generate a Perl file by copying and pasting the code below into a .pl file and run it. I hope, it will help as it helped me :-)
Code:
#!/usr/bin/perl
use strict;
use warnings;
use File::Basename;
########################
sub suffix_remover($);
################		###########################		############################		###########################
=cut
Should you decide to make this script publicly available, copy the suffix_remover() subroutine into here.
This script takes as input two fastq files with the forward and reverse unmatched mates for a paired-end read data set.
Input: provide on the command line, the paired-end read files (forward and reverse) containing reads that need to be matched. 
Output: the forward and reverse mate files will be matched and written to READNAME_sorted.fastq, where
'READNAME' is simply the name of the file provided. The READNAME_singletons.fastq files contain the singleton reads
(reads with no matching mates) for both the input files. 

by: Armin PEYMANN
=cut
################		###########################		############################		###########################
sub sequence_separator ($);
sub hash_maker_using_fastq_sequences (\@);
############################################
my $path_read1 = shift @ARGV;
my $path_read2 = shift @ARGV;
my $qx_output_read1 = qx(wc -l $path_read1);
my ($number_of_lines_read1) = split(/\s/, $qx_output_read1);
my $qx_output_read2 = qx(wc -l $path_read2);
my ($number_of_lines_read2) = split(/\s/, $qx_output_read2);
if ($number_of_lines_read1 % 4 != 0 || $number_of_lines_read2 % 4 != 0){
die "The number of lines in one of the files is not dividable by 4.\nFor each sequence in your fastq files, you must have"
		. " the following lines:\n" 
		. "\@HEADER\nSEQUENCE\n+QUALITY-HEADER\nQUALITY VALUES\n"
		. "Also make sure there is no empty line at the end of your file.\n"; 
} 


 my @read1 = @{sequence_separator($path_read1)};
 my %read1 = %{hash_maker_using_fastq_sequences(@read1)};
 undef @read1;
 my @read2 = @{sequence_separator($path_read2)};
 my %read2 = %{hash_maker_using_fastq_sequences(@read2)};
 undef @read2;
 my $dir_read1 = dirname($path_read1);
 my $read1_name_without_suffix = suffix_remover($path_read1);
 my $path_out_read1 = $dir_read1 . "/" . $read1_name_without_suffix . "_sorted.fastq";
 open(FH_OUT1, ">$path_out_read1");
 my $dir_read2 = dirname($path_read2);
 my $read2_name_without_suffix = suffix_remover($path_read2);
 my $path_out_read2 = $dir_read2 . "/" . $read2_name_without_suffix . "_sorted.fastq";
 open(FH_OUT2, ">$path_out_read2");	
 my $path_out_read1_singletons = $dir_read1 . "/" . $read1_name_without_suffix . "_singletons.fastq";
 open(FH_OUT3, ">$path_out_read1_singletons");	
 foreach my $key (sort keys %read1){	
if (exists $read2{$key}){	
	print FH_OUT1 $read1{$key};
	print FH_OUT2 $read2{$key};
}else{
	print FH_OUT3 $read1{$key}; 
}
 }
 close FH_OUT1;
 close FH_OUT2;
 close FH_OUT3;
 print "sorted reads were written to:\n";
 print "Check out $path_out_read1" . "\n";
 print "Check out $path_out_read2". "\n" . "\n";
 my $path_out_read2_singletons = $dir_read2 . "/" . $read2_name_without_suffix . "_singletons.fastq";
 open(FH_OUT4, ">$path_out_read2_singletons");
 foreach my $key (sort keys %read2){
 	unless (exists $read1{$key}){
 		print FH_OUT4 $read2{$key};
 	}
 }
 close FH_OUT4;
 print "single reads with no mate were written to:\n";
 print "Check out $path_out_read1_singletons" . "\n";
 print "Check out $path_out_read2_singletons" . "\n";
 
 sub hash_maker_using_fastq_sequences (\@){
 	my $array_ref = shift;
 	my @read = @{$array_ref};
 	my %read;
	foreach my $sequence (@read){
		my $copy_sequence = $sequence;
		my ($first_header) = split(/\n|\r/, $copy_sequence);
		$first_header =~ /(.+)[# ]/;	# one single space or '#' should cover most of fastq files.
		my $pair_id = $1;
		$read{$pair_id} = $sequence;
		
		}	
	return \%read;
 } 
 	
sub sequence_separator ($){
	my $path_read = shift;
my $line_counter = 0;
my $event;
my @events;
open(FH_IN, $path_read) or die "unable to open FH_IN1\n";
while(my $line = <FH_IN>){
	$line_counter++;
		
	if ($line_counter <= 4){		
		$event .= $line;
	}
	if ($line_counter == 4){
	push(@events, $event);
	undef $event;
	$line_counter = 0	
	}
  }
  close FH_IN;
  return \@events;
} 

sub suffix_remover($){							#get rid of the .<format> in the file name. (If there are more 
	my $pathOfDesiredFile = shift;					#than one dot in the file name, the shortest part of the file name will be taken!) 
	$pathOfDesiredFile =~ /([^\/]+)(?!\/)$/;
	my $desiredFileName = $1 if $1;
	$desiredFileName =~ s/(\..+)(?!\.)// if $desiredFileName =~ /\./ ; 														
	return $desiredFileName;	
}

Last edited by Brian Bushnell; 07-30-2015 at 08:52 AM. Reason: Wrapped in [code]
Armin_P is offline   Reply With Quote
Old 07-31-2015, 11:14 AM   #56
SES
Senior Member
 
Location: Vancouver, BC

Join Date: Mar 2010
Posts: 275
Default

Quote:
Originally Posted by Armin_P View Post
Since this has generally remained a common problem
On the contrary, I think this is a solved problem for the most part. There are at least two general solutions (Pairfq and bbtools) discussed above that work on multi-line fasta/q, compresed files, and handle different Illumina variants. The Pairfq option doesn't require any install. It might seem like this is a common issue, but in my opinion, that is because people find this thread and try the Python script posted early on and they keep running into the same issues. It appears that script is not maintained so no one can get it to work.

There's no harm in rolling your own solution, but people should know (if they haven't followed the whole discussion) that there are working (and tested, documented) solutions.
SES is offline   Reply With Quote
Old 08-21-2020, 04:00 AM   #57
jessieHOU
Junior Member
 
Location: Hong Kong

Join Date: Aug 2020
Posts: 1
Default

How to download this script?

Last edited by jessieHOU; 08-21-2020 at 04:09 AM.
jessieHOU is offline   Reply With Quote
Old 08-21-2020, 06:08 AM   #58
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,077
Default

Download BBMap suite. "repair.sh" is the program that will sync paired-end data files. It is part of BBtools ahd has a guide available.
GenoMax is offline   Reply With Quote
Reply

Tags
fastx-toolkit, paired-end reads

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 01:32 AM.


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