1
$\begingroup$

In f1 below I am trying to split $1 based on each line and create a bed file. If the line is a snp then the : is spilt and the text is $1 the last digit is $2 -1 and $3 and the letter to the left of the > is $4 and the letter to the right is $5. (lines 1 and 2)

If the line has dup then the : is spilt and the text is $1 the last digit is $2 and $3 and the $4 is - and the letter to the right is $5. (line 4)

If the line has del with nothing else then the : is spilt and the text is $1 the last digit is $2 and $3 and the letter to the left of the > is $4 and $5 is -. (line 5)

I am not sure how to format for line 3 (the complex indel)... maybe split $1 on the : and the text is $1 the last digit is $2 and $3 and the letters to the left of the > is $4 and the letters to the right is $5. Thank you :)

f1

chr7:140453145A>T
chr7:140453136A>T 
chr7:140453135_140453136delCAinsTT
chr20:31022287dupA
chr19:13054614delG

desired

chr7    140453144   140453145   A   T 
chr7    140453135   140453136   A   T 
chr7    140453135   140453136   TT  CA
chr20   31022287    31022287    -   A
chr19   13054614    13054614    G   -

awk

awk 'BEGIN{OFS="\t"} {sub(/[^0-9]+$/, "", $1)}' f1

awk 'BEGIN { FS = "[ -]"; OFS="\t" } NF==3 { print $1, $2 - 1, $2, $3 } NF==4 { print $1, $2, $3, $4 }' f1
$\endgroup$
  • 2
    $\begingroup$ Do you have to do this in awk for some reason? Sure it's possible, but there are enough cases that'd it'd make more sense to use python that's easier to write longer scripts in. $\endgroup$ – Devon Ryan Sep 23 at 12:57
  • $\begingroup$ python is fine I just don't know it as well.... only the basics but I am interested to see what it would look like as I can learn from it. Thank you :). $\endgroup$ – justaguy Sep 23 at 13:21
4
$\begingroup$

The following Python script seems to do the job.

#!/usr/bin/env python3

from __future__ import print_function
import re
import sys

for line in sys.stdin:
    if '_' in line:
        match = re.search(r'(\S+):(\d+)_(\d+)del([ACGT]+)ins([ACGT]+)', line)
        print(match.group(1), match.group(2), match.group(3), match.group(5), match.group(4), sep='\t')
    elif 'dup' in line or 'del' in line:
        match = re.search(r'(\S+):(\d+)(del|dup)([ACGT]+)', line)
        if match.group(3) == 'del':
            print(match.group(1), match.group(2), match.group(2), match.group(4), '-', sep='\t')
        else:
            print(match.group(1), match.group(2), match.group(2), '-', match.group(4), sep='\t')
    else:
        match = re.search(r'(\S+):(\d+)([ACGT]+)>([ACGT]+)', line)
        coord = int(match.group(2))
        print(match.group(1), coord - 1, coord, match.group(3), match.group(4), sep='\t')

Invoked like so on the command line.

[standage@lappy ~] $ ./transform < f1
chr7 140453144 140453145 A T
chr7 140453135 140453136 A T
chr7    140453135   140453136   TT  CA
chr20   31022287    31022287    -   A
chr19   13054614    13054614    G   -
$\endgroup$
  • $\begingroup$ Thank you very much: print(match.group(1), match.group(2), match.group(3), match.group(5), match.group(4), sep='\t') ^ SyntaxError: invalid syntax ~$ python --version Python 2.7.9 (not sure if this helps, sorry new to python. :) $\endgroup$ – justaguy Sep 23 at 14:43
  • 1
    $\begingroup$ It looks like you're running Python 2, which unfortunately is the default version on many operating systems even though it's at its end-of-life. I'll update the script so that it's compatible with Python 2. $\endgroup$ – Daniel Standage Sep 23 at 14:51
  • 1
    $\begingroup$ "print" was a statement in Python 2 but is a function in Python 3 hence the syntax error. $\endgroup$ – haci Sep 23 at 14:51
3
$\begingroup$

The Python script in my first response is heavy on regex matching, which is pretty clunky in Python. I like Python much better than Perl overall, but a throwaway script like this will be clearer and more concise in Perl.

#!/usr/bin/env perl
use strict;

while(<STDIN>) {
    if (m/(\S+):(\d+)_(\d+)del([ACGT]+)ins([ACGT]+)/) {
        print("$1\t$2\t$3\t$5\t$4\n");
    }
    elsif (m/(\S+):(\d+)(del|dup)([ACGT]+)/) {
        if ($3 eq "del") {
            print("$1\t$2\t$2\t$4\t-\n");
        }
        else {
            print("$1\t$2\t$2\t-\t$4\n");
        }
    }
    elsif (m/(\S+):(\d+)([ACGT]+)>([ACGT]+)/) {
        my $coord2 = $2;
        my $coord1 = $2 - 1;
        print("$1\t$coord1\t$coord2\t$3\t$4\n")
    }
}

It's been several years since I wrote Perl on a regular basis, so there are probably ways to make the script even more clear and concise. (Although I have a love/hate relationship with the default/hidden variables, which can make a script more concise but also harder for a newcomer to understand.)

Invoked like so on the command line.

[standage@lappy ~] $ ./transform < f1
chr7 140453144 140453145 A T
chr7 140453135 140453136 A T
chr7    140453135   140453136   TT  CA
chr20   31022287    31022287    -   A
chr19   13054614    13054614    G   -
$\endgroup$
3
$\begingroup$

Quick Perl solution (following on from Daniel's answer saying that this would be "clearer and more concise in Perl"). I was a bit confused with your BED output because you've treated the start/end locations of a single base deletion and insertion different to a SNP. In my example here, the start and end locations for those cases are the same to distinguish from the two-base change for TT->CA [this can be changed]:

cat 10454.txt | \
perl -ne 'if(/^(.*?):([0-9]+)(_([0-9]+))?(del|dup)?([ACGT]+)((ins|>)([ACGT]+))?/){
   print(join("\t",
         $1,
         $2,
         $4?$4:($2),
         ($8 eq "ins")?$9:($5 eq "dup")?"-":$6,
         (($5 eq "dup") || ($8 eq "ins"))?$6:($9)?$9:"-")."\n")
}';
chr7    140453145       140453145       A       T
chr7    140453136       140453136       A       T
chr7    140453135       140453136       TT      CA
chr20   31022287        31022287        -       A
chr19   13054614        13054614        G       -

Regular Expression

This uses enclosing groups to describe the important bits of the sequence. The groups are numbered from 1 upwards, with each number corresponding to a left bracket (in order), capturing the entirety of the sub-expression until the paired right bracket. I'll describe the groups, which are referenced in the subsequent part of the script as $<group_number>:

  1. Any number of characters, up to (but not including) the first ':'
  2. Any number of the characters '0' to '9' [i.e. numbers]
  3. Optionally, a "_", followed by any number of the characters '0' to '9'
  4. The numbers from Group [3], i.e. excluding the "_"
  5. Optionally, a "dup" or a "del"
  6. Any number of characters selected from "A", "C", "G", "T" (i.e. bases)
  7. Optionally, an "ins" or ">", followed by any number of bases
  8. The "ins" or ">" from Group [7]
  9. The bases from Group [7]

Print Code

The print statement prints out a tab-separated list, finishing off with a line break. I'll describe each list element:

  1. Group [1]
  2. Group [2]
  3. If Group [4] exists, then Group [4], otherwise Group [2]
  4. If Group [8] is "ins", then Group [9]; [otherwise]
    If Group [5] is "dup", then "-", otherwise Group [6]
  5. If Group [5] is "dup", or group [8] is "ins", then Group [6][otherwise]
    If Group [9] exists, then Group [9], otherwise "-"
$\endgroup$
  • 1
    $\begingroup$ +1, although I'm tempted to quote Jurassic Park here: "You were so preoccupied with whether or not you could, you didn't stop to think if you should." 😀 $\endgroup$ – Daniel Standage Sep 24 at 16:17
2
$\begingroup$

Here is an R solution, would probably be slower than the suggested Python solution though:

bed <- readLines("bed.txt")

res <- data.frame()

for(i in bed) {

  temp <- strsplit(i, ":") %>% unlist()

  # case N>N

  if(grepl(">", temp[[2]])) {

    begin <- str_extract(temp[[2]], pattern = "^[:digit:]+") %>% as.integer() - 1
    end <- begin + 1
    ref <- str_extract(temp[[2]], pattern = "[A-Z]")
    alt <- str_extract(temp[[2]], pattern = "[A-Z]$")

    temp_res <- data.frame(temp[[1]], begin, end, ref, alt)
    res <- rbind(res, temp_res)
  }

  # case del_ins

  if(grepl("ins", temp[[2]])) {

    begin <- str_extract(temp[[2]], pattern = "^[:digit:]+") %>% as.integer()
    end <- str_extract(temp[[2]], pattern = "_[:digit:]+")
    end <- substring(end, 2, nchar(end)) %>% as.integer()
    ref <- str_extract(temp[[2]], pattern = "[A-Z]+$")
    alt <- str_extract(temp[[2]], pattern = "[A-Z]+")

    temp_res <- data.frame(temp[[1]], begin, end, ref, alt)
    res <- rbind(res, temp_res)
  }


  # case dup

  if(grepl("dup", temp[[2]])) {

    begin <- str_extract(temp[[2]], pattern = "^[:digit:]+") %>% as.integer()
    end <- begin
    ref <- "-"
    alt <- str_extract(temp[[2]], pattern = "[A-Z]$")

    temp_res <- data.frame(temp[[1]], begin, end, ref, alt)
    res <- rbind(res, temp_res)
  }

  # case del only

  if(grepl("del", temp[[2]]) & !grepl("ins", temp[[2]])) {

    begin <- str_extract(temp[[2]], pattern = "^[:digit:]+") %>% as.integer()
    end <- begin
    ref <- str_extract(temp[[2]], pattern = "[A-Z]$")
    alt <- "-"

    temp_res <- data.frame(temp[[1]], begin, end, ref, alt)
    res <- rbind(res, temp_res)
  }

}

print(res)

  temp..1..     begin       end ref alt
1      chr7 140453144 140453145   A   T
2      chr7 140453135 140453136   A   T
3      chr7 140453135 140453136  TT  CA
4     chr20  31022287  31022287   -   A
5     chr19  13054614  13054614   G   -
$\endgroup$
  • $\begingroup$ Thank you all very much :). $\endgroup$ – justaguy Sep 23 at 20:33

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.