Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • #46
    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.

    Comment


    • #47
      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.

      Comment


      • #48
        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.

        Comment


        • #49
          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

          Comment


          • #50
            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

            Comment


            • #51
              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, 11:13 AM. Reason: fix typo

              Comment


              • #52
                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,

                Comment


                • #53
                  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.


                  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!

                  Comment


                  • #54
                    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".

                    Comment


                    • #55
                      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, 08:52 AM. Reason: Wrapped in [code]

                      Comment


                      • #56
                        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.

                        Comment


                        • #57
                          How to download this script?
                          Last edited by jessieHOU; 08-21-2020, 04:09 AM.

                          Comment


                          • #58
                            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.

                            Comment


                            • #59
                              I use BBMap v39.01 with this command:

                              repair.sh -Xmx80g t=16 in1=SG006_1000_1.fastq.gz in2=SG006_1000_1.fastq.gz out1=SG006_1000_fixed_1.fastq.gz out2=SG006_1000_1_fixed.fastq.gz outs=singletons.fq repair tossbrokenreads=t

                              and it failed with:

                              Set threads to 16
                              Set INTERLEAVED to false
                              Exception in thread "main" java.lang.AssertionError
                              at jgi.SplitPairsAndSingles.<init>(SplitPairsAndSingles.java:207)
                              at jgi.SplitPairsAndSingles.main(SplitPairsAndSingles.java:37)


                              what's wrong ?

                              Comment


                              • #60
                                mslider your error is likely an issue with the input or the processing in BBMap. One obvious issue in your command is that both `in1` and `in2` are set to the same input file (`SG006_1000_1.fastq.gz`). Typically, paired-end sequences are provided in two separate files for forward and reverse reads. Make sure `in1` and `in2` point to the correct, distinct input files.

                                Comment

                                Latest Articles

                                Collapse

                                • seqadmin
                                  Strategies for Sequencing Challenging Samples
                                  by seqadmin


                                  Despite advancements in sequencing platforms and related sample preparation technologies, certain sample types continue to present significant challenges that can compromise sequencing results. Pedro Echave, Senior Manager of the Global Business Segment at Revvity, explained that the success of a sequencing experiment ultimately depends on the amount and integrity of the nucleic acid template (RNA or DNA) obtained from a sample. “The better the quality of the nucleic acid isolated...
                                  03-22-2024, 06:39 AM
                                • seqadmin
                                  Techniques and Challenges in Conservation Genomics
                                  by seqadmin



                                  The field of conservation genomics centers on applying genomics technologies in support of conservation efforts and the preservation of biodiversity. This article features interviews with two researchers who showcase their innovative work and highlight the current state and future of conservation genomics.

                                  Avian Conservation
                                  Matthew DeSaix, a recent doctoral graduate from Kristen Ruegg’s lab at The University of Colorado, shared that most of his research...
                                  03-08-2024, 10:41 AM

                                ad_right_rmr

                                Collapse

                                News

                                Collapse

                                Topics Statistics Last Post
                                Started by seqadmin, Yesterday, 06:37 PM
                                0 responses
                                7 views
                                0 likes
                                Last Post seqadmin  
                                Started by seqadmin, Yesterday, 06:07 PM
                                0 responses
                                7 views
                                0 likes
                                Last Post seqadmin  
                                Started by seqadmin, 03-22-2024, 10:03 AM
                                0 responses
                                49 views
                                0 likes
                                Last Post seqadmin  
                                Started by seqadmin, 03-21-2024, 07:32 AM
                                0 responses
                                66 views
                                0 likes
                                Last Post seqadmin  
                                Working...
                                X