3

I have a file containing the following informations:

    gene            3025..3855
                     /gene="Sp34_10000100"
                     /ID="Sp34_10000100"
     CDS             join(3025..3106,3722..3855)
                     /gene="Sp34_10000100"
                     /codon_start=1
                     /ID="Sp34_10000100.t1.cds1,Sp34_10000100.t1.cds2"
     mRNA            3025..3855
                     /ID="Sp34_10000100.t1"
                     /gene="Sp34_10000100"
     gene            12640..13470
                     /gene="Sp34_10000200"
                     /ID="Sp34_10000200"
     CDS             join(12640..12721,13337..13470)
                     /gene="Sp34_10000200"
                     /codon_start=1
                     /ID="Sp34_10000200.t1.cds1,Sp34_10000200.t1.cds2"
     mRNA            12640..13470
                     /ID="Sp34_10000200.t1"
                     /gene="Sp34_10000200"
     gene            15959..20678
                     /gene="Sp34_10000300"
                     /ID="Sp34_10000300"
     CDS             join(15959..16080,16268..16367,18913..19116,20469..20524,20582..20678)
                     /gene="Sp34_10000300"
                     /codon_start=1
                     /ID="Sp34_10000300.t1.cds1,Sp34_10000300.t1.cds2,Sp34_10000300.t1.cds3,Sp34_10000300.t1.cds4,Sp34_10000300.t1.cds5"
     mRNA            15959..20678
                     /ID="Sp34_10000300.t1"
                     /gene="Sp34_10000300"
     gene            22255..23085
                     /gene="Sp34_10000400"
                     /ID="Sp34_10000400"

I want to delete all the gene sections but CDS and mRNA information should be there. The output should be like this:

CDS             join(3025..3106,3722..3855)
                     /gene="Sp34_10000100"
                     /codon_start=1
                     /ID="Sp34_10000100.t1.cds1,Sp34_10000100.t1.cds2"
     mRNA            3025..3855
                     /ID="Sp34_10000100.t1"
                     /gene="Sp34_10000100"
     CDS             join(12640..12721,13337..13470)
                     /gene="Sp34_10000200"
                     /codon_start=1
                     /ID="Sp34_10000200.t1.cds1,Sp34_10000200.t1.cds2"
     mRNA            12640..13470
                     /ID="Sp34_10000200.t1"
                     /gene="Sp34_10000200"
     CDS             join(15959..16080,16268..16367,18913..19116,20469..20524,20582..20678)
                     /gene="Sp34_10000300"
                     /codon_start=1
                     /ID="Sp34_10000300.t1.cds1,Sp34_10000300.t1.cds2,Sp34_10000300.t1.cds3,Sp34_10000300.t1.cds4,Sp34_10000300.t1.cds5"
     mRNA            15959..20678
                     /ID="Sp34_10000300.t1"
                     /gene="Sp34_10000300"

Please give me any suggestion how to do this.

  • It seems you have a couple of text processing tasks to do with those files. I suggest to not just ask questions and try the answers, but understand them and learn to help yourself. Then you only need to ask if you don't succeed with your own solution. By posting your own attempts we can help you best. Thanks. – Philippos Apr 20 '17 at 8:52
  • Ya, you are right. But the thing is I am new with this kind of thigns. I am also learning by myself. But I want to see how people solve this. it will help me to learn as well I think. Thanks for your valuable suggestion – Masum Billah Apr 20 '17 at 8:56
  • 1
    If all gene sections end in a /ID=" line, sed '/^ *gene/,/\/ID="/d' might be enough. – Stéphane Chazelas Apr 20 '17 at 11:11
4

awk is usually easier to read and understand:

Here is a simple program that writes by default, and toggle a "wewrite" to "0" (= off, we will not write) when it sees a line where the first word is "gene", and put it back on when he sees a line where the first word is "CDS" or "mRNA" :

awk '
  BEGIN                               { weprint=1 }

  ( $1 == "gene" )                    { weprint=0 }
  ( $1 == "CDS" ) || ( $1 == "mRNA" ) { weprint=1 }
  ( weprint == 1)                     { print $0 ;}

  '  file_to_read

BEGIN is done before any lines are read.

The other ( test ) { action if test successful } are parsed for each line of input (... unless an action contains next, which then would ignore the rest of those and instead would go fetch the next line of input)

This will only print sections "CDS" and "mRNA" and not "gene"

This could be "golfed" (for example, the default action for a successfull 'test' is to print $0, so you could have just ( weprint == 1) as the last line, but it would be less clear to grasp, imo...)

  • @masum-billah : does this work for you? I just wrote it and didn't test it (can't test at the moment, I can just write from memory) – Olivier Dulac Apr 20 '17 at 12:40
  • Yes, it worked very nicely and was pretty easy to understand by me. I am laughing now why didn"t I get the idea to use this kind of thing. – Masum Billah Apr 21 '17 at 1:34
3
sed -e '
   /^ *gene /!b   # print non-gene block begin lines
   :a  
   $d; N          # do-while loop accumulates lines for gene block
   s/\n *\///;ta
   D              # clip the gene block
' yourfile

You need to realize that the sed model is to read a file on a per-line basis, and sed command in the -e section is applied in sequence on the line as it get's transformed unless there are branching instructions involved. And a basic syntax of sed is address command where command can be any valid sed command and address can be either of these: linenum, $ (= last line), regex, range of addresses, and finally nothing meaning this gets applied to ALL lines. Note that lines are stored in a register called the pattern space.

So with that basic stuff out of the way, we go to the actual sed -e code at hand: b => branch to the end of sed code and print the pattern space. This means we keep printing any line that does NOT (the ! after the address pattern) have the string gene as it's first field.

When we finally hit the gene in the first field line, we set up a do-while loop (:a sets a mark to be jumped to) to keep accumulating the lines into the pattern space register (N appends the next line; s command removes \n *\/, which is the line break, followed by spaces and a /) till the time either of the 2 conditions are not met, viz., either we hit the eof => we delete it ($d => delete pattern space if we are at the last line) since this a gene block which appeared towards the eof and must go.

OR we hit the beginning of next block: if s could find and remove the said pattern, the t will jump to :a, otherwise (a new block, so the pattern was not found), we continue. Now the pattern space holds the whole of the gene block and the first line of the next block. We promptly delete the gene block and with the beginning of the next block we go to the top of sed code (that's what the D command does).

  • Sorry can you please bit explain it to me. I am a new learner – Masum Billah Apr 20 '17 at 8:58
  • Thanks Rakesh Sharma it worked well and thanks again for making the commands so clear to me – Masum Billah Apr 21 '17 at 1:19
1

I can't resist to give a perl answer when we have sed and awk answer!

# make perl complain when it should
use strict;
use warnings;

# declare variable
my $section;

# run through every line
while (<>) {
  # set the current section to 'gene', 'CDS' or 'mRNA' when it matches
  $section = $1 if /^\h*(gene|CDS|mRNA)/;

  # print if the current section is not 'gene'
  print if $section ne 'gene';
}
  • a welcomed addition, the code is clean and readable too, you get +1 from me – Olivier Dulac Apr 21 '17 at 10:41
  • @OlivierDulac Glad you like it! – Alex Vong Apr 23 '17 at 6:18

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.