次の9つの変数を含むtsvファイルがあります。
> seqnames start endwidth strand metadata X.10logMacsq annotation distanceToTSS
メタデータ列には分析する情報が含まれていますが、最初に項目を分割してヘッダーを持つ独自の列に入れる必要があります。メタデータは次のとおりです(最初の行)。
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;
この列には(1行あたり)合計27の項目がありますが(ここには表示されていません)、最初にその項目をすべてその列に記録してから、不要な項目を削除する必要があると思いました。説明的な列見出しがある場合は、名前を削除することもできます(例:ID = SRX SRXなど)。
サンプルファイルの入力(最初の行)
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
この問題を解決するのに役立つ、またはアドバイスをすることができる人はいますか?私はBashに初めて触れていて、これらのコマンドに慣れていません。
これまでにいくつかのファイルをクリーンアップしました。
cut --complement -f 9-14 hisHMECanno.tsv | sed 's/%20/ /g' > hisHMECannoFilt.tsv
(元のファイルには不要な列がいくつかあり、ちょうど削除しました。)
その後、awkを使用して項目をタブ区切りの列に分割しようとしましたが、役に立ちませんでした。
答え1
次のPerlスクリプトは次を使用します。テキスト::CSVモジュールはTSVファイルを読み取り、正しい形式のTSVデータを出力します。
必要に応じて自動的にフィールドを引用し、Text::CSV
設定を使用してundef_str
未定義のメタデータフィールドを引用符付きの空の文字列として出力します(またはで""
印刷する方法のコメント付きの例を含む)。N/A
--
これらの3行のうち、最大1行はコメントアウトを削除し、残りの行は削除またはコメントアウトする必要があります。このフィールドを空白のままにするには、3行すべてを削除/コメントアウトしてください。
未定義のフィールドに何かを追加することをお勧めします。これは、2つ以上のタブ(空のフィールドなど)を単一のタブと同じように扱うことができる他のツールを使用してこのスクリプトの出力を後処理するのが簡単になるためです。スペース数の制限ではなく、単一のタブに明示的に設定しない限り、デフォルトでこれを行いますawk
。perl
Text::CSV
Debianおよび関連ディストリビューション用にlibtext-csv-perl
(純粋なPerlバージョン)およびlibtext-csv-xs-perl
(より高速にコンパイルされたCモジュール)としてパッケージ化されています。を使用してくださいapt install libtext-csv-perl
。他のディストリビューションでもそれをパッケージ化できます。それ以外の場合を使用してください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;
};
たとえばprocess-tsv.pl
、これを実行可能にしてchmod +x process-tsv.pl
実行するときにファイル名引数を指定します。例えば
$ ./process-tsv.pl filename.tsv
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
もちろん、出力をシェルのファイルにリダイレクトできます。
./process-tsv.pl input.tsv > output.tsv
答え2
すべてのUnixシステムのシェルでawkを使用すると、おそらくこれを実行したいと思います。ただし、推測をテストするために使用できるサンプル入力/出力はありません。
$ 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
上記は以下の入力ファイルを使用して実行されました。
$ 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
その中のすべてのスペースはタブ文字です。
column
以下を使用してテーブルに可視化して、値がラベル(列ヘッダー文字列)とどのようにソートされるかを確認できます。
入力する:
$ 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
出力:
$ 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
答え3
使用ミラー
$ 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
これら2つのput
コマンドを組み合わせることもできますが、「データのクリーンアップ」ステップと「フィールド分割」ステップに分けることがより明確になります。
より明確なフィールド分割のために、出力形式をCSVに変更します。
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
答え4
これは、正規表現とリストの理解と組み合わせたPython辞書とリストのデータ構造を使用して実行できます。
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
出力:
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