Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • BasesBreaker
    replied
    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.

    Leave a comment:


  • mslider
    replied
    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 ?

    Leave a comment:


  • GenoMax
    replied
    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.

    Leave a comment:


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

    Leave a comment:


  • SES
    replied
    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.

    Leave a comment:


  • Armin_P
    replied
    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]

    Leave a comment:


  • Brian Bushnell
    replied
    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".

    Leave a comment:


  • SES
    replied
    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!

    Leave a comment:


  • luqmanaslam
    replied
    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,

    Leave a comment:


  • SES
    replied
    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

    Leave a comment:


  • luqmanaslam
    replied
    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

    Leave a comment:


  • Brian Bushnell
    replied
    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

    Leave a comment:


  • SES
    replied
    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.

    Leave a comment:


  • Brian Bushnell
    replied
    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.

    Leave a comment:


  • SES
    replied
    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.

    Leave a comment:

Latest Articles

Collapse

  • seqadmin
    Best Practices for Single-Cell Sequencing Analysis
    by seqadmin



    While isolating and preparing single cells for sequencing was historically the bottleneck, recent technological advancements have shifted the challenge to data analysis. This highlights the rapidly evolving nature of single-cell sequencing. The inherent complexity of single-cell analysis has intensified with the surge in data volume and the incorporation of diverse and more complex datasets. This article explores the challenges in analysis, examines common pitfalls, offers...
    06-06-2024, 07:15 AM
  • seqadmin
    Latest Developments in Precision Medicine
    by seqadmin



    Technological advances have led to drastic improvements in the field of precision medicine, enabling more personalized approaches to treatment. This article explores four leading groups that are overcoming many of the challenges of genomic profiling and precision medicine through their innovative platforms and technologies.

    Somatic Genomics
    “We have such a tremendous amount of genetic diversity that exists within each of us, and not just between us as individuals,”...
    05-24-2024, 01:16 PM

ad_right_rmr

Collapse

News

Collapse

Topics Statistics Last Post
Started by seqadmin, Today, 06:54 AM
0 responses
10 views
0 likes
Last Post seqadmin  
Started by seqadmin, 06-14-2024, 07:24 AM
0 responses
16 views
0 likes
Last Post seqadmin  
Started by seqadmin, 06-13-2024, 08:58 AM
0 responses
15 views
0 likes
Last Post seqadmin  
Started by seqadmin, 06-12-2024, 02:20 PM
0 responses
17 views
0 likes
Last Post seqadmin  
Working...
X