Take the 2-minute tour ×
Biology Stack Exchange is a question and answer site for biology researchers, academics, and students. It's 100% free, no registration required.

I am trying to check if a fastq file has single or paired end reads. How can I achieve this with an error-proof method?

I checked wikipedia and MAQ but I want to know if is there a reliable document that describes all possible variants in sequence ID to check for single/paired end reads.

I am searching also for a library, better in Python, to achieve this.

Thanks

share|improve this question
    
You might be better off asking this question on biostars.org –  jarlemag Mar 5 at 17:34
    
See for example this question: biostars.org/p/45961 –  jarlemag Mar 5 at 17:49
    
Thanks @jarlemag. I did some research earlier and I noticed that answer. However there are two problems: first I am not fluent in Perl and it is hard to me to deciphers the content of the accepted answer (not to mention that not being able to check the code could lead to run into bugs that you cannot predict). Second, I read the answer and wikipedia page for FASTQ format but I am searching for the real reference. So I decided to ask a new question not to have a script (related but not main question) but the idea behind that, that could possibly cover all possible variants in this scenario. –  nsl Mar 5 at 23:40
1  
I see. I found some info on the read identifiers, but held off from answering since it seemed you wanted a complete script. Have you checked out the SeqIO interface in BioPython? –  jarlemag Mar 5 at 23:47
    
Not yet, sorry.. I am quite a noob in BioPython.. now I'll try it. Thanks –  nsl Mar 5 at 23:49

1 Answer 1

By now I got some interesting answers in this question on Biostars

Basically what I did is the following:

  • First of all, I checked if Sequence Id contains paired end notation. As described in this wikipedia page, for Illumina reads there are two possible notation for single/paired end reads:

    @HWUSI-EAS100R:6:73:941:1973#0/1

    If the last number is /2 in some reads then the reads are paired end; otherwise they can be single end.

    The second notation is:

    @EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG

    If the first number in the second group is 2 in some reads then the reads are paired end; otherwise they can be single end;

  • Then I checked for multiple files. If a sample has two fastq files it is likely that the reads are paired end. It is to note however that with a single file it is not possible to exclude that paired end reads can be interleaved in one single file, even if it is not common (in my opinion);

  • The most general method is to cross check each single reads with the whole set of reads. If the first part of the sequence ID (in this case the field starting from @ and ending before the # - in the first notation - or the whitespace - in the second notation) is unique among all the reads (for each read) it is likely that the reads are single reads, otherwise - if can be found a duplicate for each reads - the reads are paired end. In this case, on *nix systems it can be achieved with the following command (thanks to the biostars answers):

    grep --no-filename @HWUSI-EAS100R:6:73:941:1973 *.fastq | cut -d' ' -f1 | sort | uniq -c | sort -rgk 1,1 | head

    If the result shows in the first lines a result like this:

    1 read1_ID

    1 read2_ID

    ...

    It is likely it is single end. Otherwise:

    2 read1_ID

    2 read2_ID

    ...

    it is paired end.

I skim read BioPython API documentation but I cannot find something useful to do it.

Suggestions and corrections are welcomed.

Thanks

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.