1

I'm trying to parse a large file where every 2 consecutive lines have the same length (the text is completely different). I've searched, and my first post here. I found a script and tried to modify it but no joy. File is a sequencing output file. I have already parsed out the sequence, and quality scores, so the file looks like this:

CCTCGNAACCCAAAAACTTTGATTTCTNATAAGGTGCCAGCGGAGTCCTAAAAGCAACATCCGCTGATCCCTGGT
AAAAA#EEEEEEEEEEEEEEEEEEEEE#EEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
CCCCANCCAAACTCCCCACCTGACAATNTCCTCCGCCCGGATCGACCCGCCGAAGCGAGTCTTGGGTCTAAA
AAAAA#EEEEEEEEEEEAEEEEEEEEE#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
ATCGTNTATGGTTGAGACTAGGACGGTNTCTGATCGTCTTCGAGCCCCCAACTTTCGTTCTTGATTAATGAAAAC
AAAAA#EEEEEEEEEEEEEEEEEEEEE#EEEEEEEEEEEEAEEEEEEAEEEAEEEEEEEEEEEEEEEEEEEEEEE
CCCACNTGGAGCTCTCGATTCCGTGGGNTGGCTCAACAAAGCAGCCACCCCGTCCTACCTATTTAAAGTTTG
AAAAA#EEEEEEEEEEEEEEEEEEEEE#EEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEE
GCATCNTTTATGGTTGAGACTAGGACGNTATCTGATCGTCTTCGAGCCCCCAACTTTCGTTCTTGATTAATGAA
6AAAA#EEEEEAAAEEEEEEAEEAEEE#EEEEEEEAEAEEEEAEEAAA/EAEEEEAEEAEEAEEAEAAEEEEEE

The problem: Somewhere there is a corrupted pair of lines such that each sequence base, does not have a corresponding score, i.e. the lengths of every pair of two lines should be equal, how can I parse out the pair that is incorrect? File is 100-million lines.

I tried this code named parser.sh:

{ curr = $0 }
(NR%2)==0 {
    currLgth = length(curr)
    prevLgth = length(prev)
    maxLgth = (currLgth > prevLgth ? currLgth : prevLgth)
    if (prevLgth==currLgth) {
        print ""
        print prevLgth
        print currLgth
        for (i=1; i<=maxLgth; i++) {
        }
    }
}
{ prev = curr }

and would run awk -f parser.sh filename but this printed out all the lines lengths even though I was using "not equal" ('==').

75
75

72
72

75
75

72
72

Am not a coder, so apologies in advance, but need help with this. Usually can find code and modify it to work, but not in this instance. -p

Fastq files have four lines for each read. Read#1 e,g, will have the following 4 lines:

@sample1
CGGCATCGTTTATGGTTGAGACTAGGACG
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEE

The first line is the sample name, the second line is the actual sequence, the third line is a '+' symbol and the fourth line is a set of ASCII "scores" for each base in the sequence. Each base has exactly one score, hence the length of line 2 must equal the length of line four. I had parsed out lines 2 and 4, looking for pairs of line with unequal length. Instead I got what looks like the pairing was lost.

Here's an example of what a FASTQ file might look like, with the question marks representing the lost or unparsed quality scores:

@sample1
CGGCATCGTTTATGGTTGAGACTAGGACG
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEE
@sample2
CCGGCTTCCGGTTCATCCCGCATCGCCAGTTC
+
@sample3
AAAA6E6/EEEEEEEE6/EE/EEAEEAA//E/
+
@sample4
ATTTCGGGGGGGGGGGGGG
+
??????????????????????????????????
@Sample5
GGTTAGCGCGCAGTTGGGCACCGTAACCCGGCTT
+
AAAAAEEEEEAEEEEEEEEEEEEEEEEEE//<EE
@sample6
CTAACCTGTCTCACGACGGTCTAAACCCAGCTCA
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEE

Here is what my (line2 + line4) parsed files looked like:

CGGCATCGTTTATGGTTGAGACTAGGACG
AAAAAEEEEEEEEEEEEEEEEEEEEEEEE
CCGGCTTCCGGTTCATCCCGCATCGCCAGTTC
AAAA6E6/EEEEEEEE6/EE/EEAEEAA//E/
ATTTCGGGGGGGGGGGGGG
GGTTAGCGCGCAGTTGGGCACCGTAACCCGGCTT
AAAAAEEEEEAEEEEEEEEEEEEEEEEEE//<EE
CTAACCTGTCTCACGACGGTCTAAACCCAGCTCA
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEE

There are consecutive two sequence lines with no quality score line between them:

ATTTCGGGGGGGGGGGGGG
GGTTAGCGCGCAGTTGGGCACCGTAACCCGGCTT

Using the code you gave me:

awk 'NR%2==0 && length($0)!=last{print "Bad pair at lines",NR-1,"and",NR}{last=length($0)}' Fastq-seq-qual-parsed.txt
Bad pair at lines 5 and 6

OR: ./new-try.awk

  • how can I parse out the pair that is incorrect - what should be the expected result for your current input fragment? – RomanPerekhrest Apr 4 '18 at 18:52
  • 3
    I would suggest that an .sh extension on a file containing awk code is confusing. – glenn jackman Apr 4 '18 at 18:55
  • Sorry for the poor filenames. – hoytpr Apr 4 '18 at 19:29
  • The FastQ file that you show has partial data for some records. @sample2 and @sample3 are missing scores completely. Is this what you're working with? – Kusalananda Apr 5 '18 at 17:47
3

I would suggest

awk '
    { first = $0; getline; second = $0 }
    length(first) != length(second) {
        print "Error at line", NR-1
        print first
        print second
    }
' file

Could use plain bash too, but it will be much slower:

nr=1
while IFS= read -r first; IFS= read -r second; do 
    if (( ${#first} != ${#second} )); then 
        printf "%s\n" "problem at line $nr" "$first" "$second"
    fi
    ((nr+=2))
done < file
1

Try:

awk 'NR%2==0 && length($0)!=last{print "Bad pair at lines",NR-1,"and",NR} {last=length($0)}' file

Example

Let's take this as a test file:

$ cat file
good123
good345
bad12
bad123
good_again
good_also1

Using our command, the unmatched pair is correctly identified:

$ awk 'NR%2==0 && length($0)!=last{print "Bad pair at lines",NR-1,"and",NR} {last=length($0)}' file
Bad pair at lines 3 and 4

How it works

  • NR%2==0 && length($0)!=last{print "Bad pair at lines",NR-1,"and",NR}

    When we are on an even numbered line, NR%2==0, we check to see if the line length is the same as that of the previous line. If it is not the same, length($0)!=last, we print a message.

  • last=length($0)

    This saves the length of the current line in variable last.

Multiple line version

For those who prefer their code spread out over multiple lines:

awk '
    NR%2==0 && length($0)!=last {
        print "Bad pair at lines",NR-1,"and",NR
    }

    {
        last=length($0)
    }' file

How to print specific lines from a file

To print, for example, line 3 from a file, we can use:

$ awk 'NR==3' file
bad12

To print a range, say all lines from 3 to 6, we can use:

$ awk 'NR>=3 && NR<=6' file
bad12
bad123
good_again
good_also1

Alternatively, we can get similar results from sed using:

$ sed -n '3p' file
bad12
$ sed -n '3,6p' file
bad12
bad123
good_again
good_also1

Using unfiltered input data

Consider this input file:

$ cat File
@sample1
CGGCATCGTTTATGGTTGAGACTAGGACG
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEE
@sample2
CCGGCTTCCGGTTCATCCCGCATCGCCAGTTC
+
@sample3
AAAA6E6/EEEEEEEE6/EE/EEAEEAA//E/
+
@sample4
ATTTCGGGGGGGGGGGGGG
+
??????????????????????????????????
@Sample5
GGTTAGCGCGCAGTTGGGCACCGTAACCCGGCTT
+
AAAAAEEEEEAEEEEEEEEEEEEEEEEEE//<EE
@sample6
CTAACCTGTCTCACGACGGTCTAAACCCAGCTCA
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
@sample7
CTAACCTGTCTCACGACGGTCTAAACCCAGCTCA
+
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEE

We can detect bad samples, meaning samples with either unequal line lengths or second lines that start with ?, as follows:

$ awk '/^\+/{next} /^@/{s=$0;n=NR;next} prev{if(/^\?/ || length(prev)!=length($0)) printf "Sample %s (line %s) is bad:\n%s\n%s\n",s,n,prev,$0;prev="";next} {prev=$0}' File
Sample @sample4 (line 11) is bad:
ATTTCGGGGGGGGGGGGGG
??????????????????????????????????
Sample @sample7 (line 23) is bad:
CTAACCTGTCTCACGACGGTCTAAACCCAGCTCA
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEE

Alternatively, if we want to ignore samples whose second ('quality') line starts with ?, then:

$ awk '/^\+/{next} /^@/{s=$0;n=NR;next} prev{if(!/^\?/ && length(prev)!=length($0)) printf "Sample %s (line %s) is bad:\n%s\n%s\n",s,n,prev,$0;prev="";next} {prev=$0}' File
Sample @sample7 (line 23) is bad:
CTAACCTGTCTCACGACGGTCTAAACCCAGCTCA
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEE
  • Thanks for the excellent annotation of your script. it also worked perfectly and printed out the lines which helped show at least one problem with the file. Both your script and glenn jackman's script ran until they hit lines 5502951 and 5502952. at that point, there were two consecutive sequence lines, with no quality lines. Everything was out of phase at that point and there were thousands of errors. – hoytpr Apr 4 '18 at 19:56
  • Now, I need to grep out that region to verify. Can I grep out a specific line? How can I print out that region of the file? – hoytpr Apr 4 '18 at 20:00
  • The top of the output file from @John1024 is: '>head newtryfirst.txt Error at line 5502951 GCCGCCCCAGCCAAACTCCCCACCTGACAATGTCCTCCGCCCGGATCGACCCGCCGAAGCGAGTCTTGGGTCTA TGAATACATCAGTGTAGCGCGCGTGCGGCCCAGAACATCTAAGGGCATCACAGACCTGTTATTGCCTCAAACTTC Error at line 5502959 AAA/AEEE/AAEEAEAE6EEEEEEEAEAEEEEEEEEEEAEEEE/EEEEEEEAEE<EE6E/EAEEEEAEEAAE<EE CCCGTCGCGTATTTAAGTCGTCTGCAAAGGATTCTACCCGCCACTCGGTGGTAATTATAATTCAAGGCGGTCCG' – hoytpr Apr 4 '18 at 20:04
  • @hoytpr To print a specific line, say line 5502951, from a file, one can use awk 'NR==5502951' file. See the new section at the end of my updated answer for more. – John1024 Apr 4 '18 at 21:43
  • I cannot reproduce the out-of-phase problem with either my code or @GlennJackman's code. Can you provide us with some minimal example that would show us this behavior? – John1024 Apr 4 '18 at 21:47
0

First make a test file where lines 5 and 6 have unequal lengths so there's something to find, (the "cccc ccc" below):

printf '%s\n' aaa aaa bbbb bbbb cccc ccc ddd ddd > foo

Abstract foo into two virtual files using bash process substitution and sed, where every char is replaced by a .:

  • the first virtual file abstracts the actual file,
  • the 2nd virtual file abstracts only the odd lines, which are then duplicated -- so that in the 2nd file each consecutive odd and even line has the same length.

...then diff those files:

diff <(sed 's/././g' foo) <(sed -n '1~2{s/././g;p;p}' foo)

Output shows line 6 doesn't match:

6c6
< ...
---
> ....

If the above output is too verbose, diff and kindred programs have many options, or can be filtered as needed. To show just the line numbers:

diff <(sed 's/././g' foo) <(sed -n '1~2{s/././g;p;p}' foo) | 
sed -n 's/c.*//p'

Output:

6

Or with more detail, i.e. numbered non-matching original file lines:

f=foo
diff <(sed 's/././g' $f) <(sed -n '1~2{s/././g;p;p}' $f) |  
sed -n 's/^\(.*\)c.*/\1/p' | grep -B 1 -wf - <(cat -n $f)

Output:

     5  cccc
     6  ccc

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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