0

I have a tsv file with 9 variables, as follows:

> seqnames  start   endwidth    strand  metadata    X.10logMacsq    annotation  distanceToTSS

The metadata column contains information I want to do some analysis on, but I first need to split the entries and put them into their own columns (with a heading). The metadata looks like this (1st row):

ID=SRX067411;Name=H3K27ac%20(@%20HMEC);Title=GSM733660:%20Bernstein%20HMEC%20H3K27ac;Cell%20group=Breast;<br>source_name=HMEC;biomaterial_provider=Lonza;lab=Broad;lab%20description=Bernstein%20-%20Broad%20Institute;datatype=ChipSeq;datatype%20description=Chromatin%20IP%20Sequencing;cell%20organism=human;cell%20description=mammary%20epithelial%20cells;cell%20karyotype=normal;cell%20lineage=ectoderm;cell%20sex=U;antibody%20antibodydescription=rabbit%20polyclonal.%20Antibody%20Target:%20H3K27ac; 

There are 27 entries for this column (per row) in total (not all shown here) but I thought I'd write them all into their own columns first, and then remove the ones I don't need afterwards. Once they have a descriptive column heading, then I could also remove their names (e.g.: ID=SRX would simply be SRX and so forth)

Sample file input (1st row)

seqnames    start   end width   strand  metadata    X.10logMacsq    annotation  geneChr geneStart   geneEnd geneLength  geneStrand  geneId  distanceToTSS
chr2    1711333 1711568 236 *   ID=SRX067411;Name=H3K27ac%20(@%20HMEC);Title=GSM733660:%20Bernstein%20HMEC%20H3K27ac;Cell%20group=Breast;<br>source_name=HMEC;biomaterial_provider=Lonza;lab=Broad;lab%20description=Bernstein%20-%20Broad%20Institute;datatype=ChipSeq;datatype%20description=Chromatin%20IP%20Sequencing; 447 Intron (uc002qxa.3/7837, intron 1 of 22)    1   1635659 1748291 112633  2   7837    36723

Can anyone help me out with this or give some advice? I am fairly new to Bash and not very comfortable with the commands yet.

So far I've just managed to clean the file up a bit with:

cut --complement -f 9-14 hisHMECanno.tsv | sed 's/%20/ /g' > hisHMECannoFilt.tsv

(the orginal files had some unneccessary columns that I just removed as well)

Then I have been trying to use awk to separate the entries into tab separated columns, but to no avail.

0

4 Answers 4

1

The following perl script uses the Text::CSV module to read the TSV file and to output properly formatted TSV data.

It automatically quotes fields if required and uses Text::CSV's undef_str setting to output undefined metadata fields as a quoted empty string "" (with commented out examples of how to print them as N/A or -- instead).

At most, only one of those 3 lines should be uncommented, the others should be deleted or commented out. If you just want those fields to be empty, delete/comment out all three of those lines.

I recommend having something in those undefined fields because it will make it easier to post-process the output of this script with other tools that might treat two or more tabs (i.e. an empty field) the same as just a single tab (e.g. both awk and perl will do that by default, unless you tell them not to by explicitly setting the field separator to a single tab, rather than the default of "any amount of whitespace").

Text::CSV is packaged for debian and related distros as libtext-csv-perl (pure perl version) and libtext-csv-xs-perl (faster compiled C module). Install both with apt install libtext-csv-perl. Other distros probably have it packaged too. Otherwise, install it with cpan.

#!/usr/bin/perl

use strict;
use Text::CSV qw(csv);

my $csv=Text::CSV->new({sep_char => "\t", quote_space => 0});

# optional: define how to print undefined fields
#$csv->undef_str ('--');
#$csv->undef_str ('N/A');
$csv->undef_str ('""');

# get header line, split into an arrayref called $cols
my $cols = $csv->getline(*ARGV);

# get first data row, extract headers & data from metadata field
my $row = $csv->getline(*ARGV);

# The following line assumes that the metadata in the FIRST data row
# contains ALL of the metadata fields in the exact order you want them
# included in the output.
#
my $md_headers = extract_metadata_headers($$row[4]);
#
# If this is not the case, then delete the extract_metadata_headers
# subroutine and define the metadata fields manually with something
# like:
#
#my $md_headers = [
#  'ID', 'Name', 'Title', 'Cell group', 'source_name',
#  'biomaterial_provider', 'lab', 'lab description', 'datatype',
#  'datatype description', 'cell organism', 'cell description',
#  'cell karyotype', 'cell lineage', 'cell sex',
#  'antibody antibodydescription'
#];
# This defines both the extra metadata headers **and** the order
# that they will be included in each output row.

# extract the data from the metadata field
my $md_data = extract_metadata($$row[4]);

# replace the metadata header in $cols aref with the md headers
splice @$cols,4,1,@$md_headers;

# replace the metadata field in $row aref with the md fields
splice @$row,4,1,@$md_data;

# print the updated header line and the first row of data
$csv->say(*STDOUT,$cols);
$csv->say(*STDOUT,$row);

# main loop: extract and print the rest of the data
while (my $row = $csv->getline(*ARGV)) {
  my $md_data = extract_metadata($$row[4]);
  splice @$row,4,1,@$md_data;

  $csv->say(*STDOUT,$row);
}

###
### subroutines
###

sub extract_metadata_headers {
  my $md = clean_metadata(shift);
  my @metadata = split /;/, $md;
  my @headers=();

  foreach (@metadata) {
    next if m/^\s*$/; # skip empty metadata
    my ($key,$val) = split /=/;
    push @headers, $key;
  };

  return \@headers;
};

sub extract_metadata {
  my $md = clean_metadata(shift);
  my @metadata = split /;/, $md;
  my %data=();

  foreach (@metadata) {
    next if m/^\s*$/; # skip empty metadata
    my ($key,$val) = split /=/;
    $data{$key} = $val;
  };

  return [@data{@$md_headers}];
};

sub clean_metadata {
    my $md = shift;
    $md =~ s/%(\d\d)/chr hex $1/eg; # decode %-encoded spaces etc.
    $md =~ s/<[^>]*>//g;            # remove HTML crap like <br>
    return $md;
};

Save it as, e.g., process-tsv.pl, make it executable with chmod +x process-tsv.pl and give it a filename argument when you run it. e.g.

$ ./process-tsv.pl filename.tsv

It will produce output like this to stdout:

$ ./process-tsv.pl input.tsv
seqnames        start   endwidth        strand  ID      Name    Title   Cell group      source_name     biomaterial_provider    lab     lab description datatype        datatype description    cell organism   cell description      cell karyotype   cell lineage    cell sex        antibody antibodydescription    X.10logMacsq    annotation      distanceToTSS
seq1    1       10      X       SRX067411       H3K27ac (@ HMEC)        GSM733660: Bernstein HMEC H3K27ac       Breast  HMEC    Lonza   Broad   Bernstein - Broad Institute     ChipSeq Chromatin IP Sequencing human   mammary epithelial cells       normal  ectoderm        U       rabbit polyclonal. Antibody Target: H3K27ac     x10     annot   dist
seq2    2       20      Y       SRX067411       H3K27ac (@ HMEC)        GSM733660: Bernstein HMEC H3K27ac       ""      ""      Lonza   Broad   Bernstein - Broad Institute     ChipSeq Chromatin IP Sequencing human   mammary epithelial cells       normal  ectoderm        U       ""      Y10     annot2  dist2

You can, of course, redirect the output to a file in your shell:

./process-tsv.pl input.tsv > output.tsv
1

Using any awk in any shell on every Unix box, this might be what you're trying to do but without sample input/output we can test with it's a guess:

$ cat tst.awk
BEGIN { FS=OFS="\t" }
{
    gsub(/%20/," ")
    gsub(/<br>/,"")
}
NR==1 {
    hdr = $0
    next
}
NR==2 {
    orig = $0
    gsub(/=[^=;]+;/,OFS,$6)
    sub(OFS"$","",$6)
    tags = $6
    $0 = hdr
    $6 = tags
    print
    $0 = orig
}
{
    gsub(/[^=;]+=/,OFS,$6)
    sub("^"OFS,"",$6)
    gsub(/;/,"",$6)
    print
}

$ awk -f tst.awk file
>       seqnames        start   endwidth        strand  ID      Name    Title   Cell group      source_name     biomaterial_provider    lab     lab description datatype        datatype description    cell organism   cell description        cell karyotype  cell lineage    cell sex        antibody antibodydescription    X.10logMacsq    annotation      distanceToTSS
>       foo     bar     etc     anon    SRX067411       H3K27ac (@ HMEC)        GSM733660: Bernstein HMEC H3K27ac       Breast  HMEC    Lonza   Broad  Bernstein - Broad Institute      ChipSeq Chromatin IP Sequencing human   mammary epithelial cells        normal  ectoderm        U       rabbit polyclonal. Antibody Target: H3K27ac     end     more    stuff

The above was run using this input file:

$ cat file
>       seqnames        start   endwidth        strand  metadata        X.10logMacsq    annotation      distanceToTSS
>       foo     bar     etc     anon    ID=SRX067411;Name=H3K27ac%20(@%20HMEC);Title=GSM733660:%20Bernstein%20HMEC%20H3K27ac;Cell%20group=Breast;<br>source_name=HMEC;biomaterial_provider=Lonza;lab=Broad;lab%20description=Bernstein%20-%20Broad%20Institute;datatype=ChipSeq;datatype%20description=Chromatin%20IP%20Sequencing;cell%20organism=human;cell%20description=mammary%20epithelial%20cells;cell%20karyotype=normal;cell%20lineage=ectoderm;cell%20sex=U;antibody%20antibodydescription=rabbit%20polyclonal.%20Antibody%20Target:%20H3K27ac;       end     more    stuff

where all spaces are tabs.

You can see how the values line up with the tags (column header strings) by using column to make it visually tabular:

Input:

$ column -s$'\t' -t file
>  seqnames  start  endwidth  strand  metadata                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           X.10logMacsq  annotation  distanceToTSS
>  foo       bar    etc       anon    ID=SRX067411;Name=H3K27ac%20(@%20HMEC);Title=GSM733660:%20Bernstein%20HMEC%20H3K27ac;Cell%20group=Breast;<br>source_name=HMEC;biomaterial_provider=Lonza;lab=Broad;lab%20description=Bernstein%20-%20Broad%20Institute;datatype=ChipSeq;datatype%20description=Chromatin%20IP%20Sequencing;cell%20organism=human;cell%20description=mammary%20epithelial%20cells;cell%20karyotype=normal;cell%20lineage=ectoderm;cell%20sex=U;antibody%20antibodydescription=rabbit%20polyclonal.%20Antibody%20Target:%20H3K27ac;  end           more        stuff

Output:

$ awk -f tst.awk file | column -s$'\t' -t
>  seqnames  start  endwidth  strand  ID         Name              Title                              Cell group  source_name  biomaterial_provider  lab    lab description              datatype  datatype description     cell organism  cell description          cell karyotype  cell lineage  cell sex  antibody antibodydescription                 X.10logMacsq  annotation  distanceToTSS
>  foo       bar    etc       anon    SRX067411  H3K27ac (@ HMEC)  GSM733660: Bernstein HMEC H3K27ac  Breast      HMEC         Lonza                 Broad  Bernstein - Broad Institute  ChipSeq   Chromatin IP Sequencing  human          mammary epithelial cells  normal          ectoderm      U         rabbit polyclonal. Antibody Target: H3K27ac  end           more        stuff
0
0

Using Miller

$ mlr --tsv put -S '
    $metadata = gsub($metadata,"%20"," "); $metadata = gsub($metadata,"<br>|;$","")
  ' then put -S '
  $* = mapsum($*,splitkvx($metadata,"=",";"))
' then cut -x -f metadata HisHMECanno.tsv
seqnames        start   end     width   strand  X.10logMacsq    annotation      geneChr geneStart       geneEnd geneLength                                                 geneStrand       geneId  distanceToTSS   ID      Name    Title   Cell group      source_name     biomaterial_provider   lab                                                 lab description  datatype        datatype description
chr2    1711333 1711568 236     *       447     Intron (uc002qxa.3/7837, intron 1 of 22)        1       1635659 1748291                                                    112633   2       7837    36723   SRX067411       H3K27ac (@ HMEC)        GSM733660: Bernstein HMEC H3K27ac       Breast HMEC                                                Lonza    Broad   Bernstein - Broad Institute     ChipSeq Chromatin IP Sequencing

The two put commands could be combined, however I think it's clearer to separate them into a "data cleanup" step and a "field splitting" step.

Changing the output format to CSV to make the field splits clearer:

mlr --itsv --ocsv put -S '
  $metadata = gsub($metadata,"%20"," "); $metadata = gsub($metadata,"<br>|;$","")
' then put -S '
  $* = mapsum($*,splitkvx($metadata,"=",";"))
' then cut -x -f metadata HisHMECanno.tsv
seqnames,start,end,width,strand,X.10logMacsq,annotation,geneChr,geneStart,geneEnd,geneLength,geneStrand,geneId,distanceToTSS,ID,Name,Title,Cell group,source_name,biomaterial_provider,lab,lab description,datatype,datatype description
chr2,1711333,1711568,236,*,447,"Intron (uc002qxa.3/7837, intron 1 of 22)",1,1635659,1748291,112633,2,7837,36723,SRX067411,H3K27ac (@ HMEC),GSM733660: Bernstein HMEC H3K27ac,Breast,HMEC,Lonza,Broad,Bernstein - Broad Institute,ChipSeq,Chromatin IP Sequencing
0

This can be done using Python dictionary and list data structures., in conjunction with regexes and list comprehension.

python3 -c 'import sys, re
ifile = sys.argv[1]
fs,rs = ofs,ors = "\t","\n"
with open(ifile) as f:
  for nr,l in enumerate(f,1):
    F = l.rstrip(rs).split(fs)
    if nr == 1:
      H = F
      idx_md = F.index("metadata")
      continue
    md_hdrs = re.findall(r"[^=;]+(?==)",F[idx_md])
    md = dict(t.split("=") for t in re.sub(r";+$","",F[idx_md]).split(";"))
    if nr == 2:
      print(*H[:idx_md], *md_hdrs, *H[idx_md+1:], sep=ofs)
    print(*F[:idx_md], *[md.get(key,"") for key in md_hdrs], *F[idx_md+1:], sep=ofs)
' file

Output:

seqnames    start   end width   strand  ID  Name    Title   Cell%20group    <br>source_name biomaterial_provider    lab lab%20description   datatype    datatype%20description  X.10logMacsq    annotation  geneChr geneStart   geneEnd geneLength  geneStrand  geneId  distanceToTSS
chr2    1711333 1711568 236 *   SRX067411   H3K27ac%20(@%20HMEC)    GSM733660:%20Bernstein%20HMEC%20H3K27ac Breast  HMEC    Lonza   Broad   Bernstein%20-%20Broad%20Institute   ChipSeq Chromatin%20IP%20Sequencing 447 Intron (uc002qxa.3/7837, intron_1_of_22)    1   1635659 1748291 112633  2   7837    36723

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .