我有一个大GTF文件,如下所示:
# ./stringtie -p 4 -G /home/humangenome_hg19/homo_gtf_file.gtf -o strAD1_as/transcripts.gtf -l strAD1 /home/software/star-2.5.2b/bin/Linux_x86_64/mapA1Aligned.sortedByCoord.out.bam
# StringTie version 1.3.2d
1 StringTie transcript 30267 31109 1000 + . gene_id "strAD1.1"; transcript_id "strAD1.1.1"; reference_id "ENST00000469289"; ref_gene_id "ENSG00000243485"; ref_gene_name "MIR1302-10"; cov "0.028725"; FPKM "0.053510"; TPM "0.109957";
1 StringTie exon 30267 30667 1000 + . gene_id "strAD1.1"; transcript_id "strAD1.1.1"; exon_number "1"; reference_id "ENST00000469289"; ref_gene_id "ENSG00000243485"; ref_gene_name "MIR1302-10"; cov "0.014218";
1 StringTie exon 30976 31109 1000 + . gene_id "strAD1.1"; transcript_id "strAD1.1.1"; exon_number "2"; reference_id "ENST00000469289"; ref_gene_id "ENSG00000243485"; ref_gene_name "MIR1302-10"; cov "0.072139";
我想要第 9 列只有gene_id
, transcript_id
,reference_id
和ref_gene_id
。它们位于第 9 列并用空格分隔(列本身是制表符分隔的)。您能帮我看看如何在 Linux 中使用简单的命令来创建这样的专栏吗?我不想使用 Excel。
答案1
理想情况下,由于数据是 GTF 格式,因此应该使用 GTF 解析器来解析它。我目前没有安装这样的解析器或解析库,因此我的解决方案仅基于您在问题中提供的数据。
要提取第 9 列:
$ cut -f 9 data.gtf
gene_id "strAD1.1"; transcript_id "strAD1.1.1"; reference_id "ENST00000469289"; ref_gene_id "ENSG00000243485"; ref_gene_name "MIR1302-10"; cov "0.028725"; FPKM "0.053510"; TPM "0.109957";
gene_id "strAD1.1"; transcript_id "strAD1.1.1"; exon_number "1"; reference_id "ENST00000469289"; ref_gene_id "ENSG00000243485"; ref_gene_name "MIR1302-10"; cov "0.014218";
gene_id "strAD1.1"; transcript_id "strAD1.1.1"; exon_number "2"; reference_id "ENST00000469289"; ref_gene_id "ENSG00000243485"; ref_gene_name "MIR1302-10"; cov "0.072139";
为了从中获得我们想要的数据,我们需要分别处理转录本和外显子,因为它们的属性在数据中具有不同的顺序。我们awk
根据当前行是否包含字符串来执行此操作并输出输入数据中的不同字段exon_number
:
$ cut -f 9 data.gtf | awk '/exon_number/ { print $2, $4, $8, $10; next } { print $2, $4, $6, $8 }'
"strAD1.1"; "strAD1.1.1"; "ENST00000469289"; "ENSG00000243485";
"strAD1.1"; "strAD1.1.1"; "ENST00000469289"; "ENSG00000243485";
"strAD1.1"; "strAD1.1.1"; "ENST00000469289"; "ENSG00000243485";
然后我们从中删除双引号和分号:
$ cut -f 9 data.gtf | awk '/exon_number/ { print $2, $4, $8, $10; next } { print $2, $4, $6, $8 }' | tr -d '";'
strAD1.1 strAD1.1.1 ENST00000469289 ENSG00000243485
strAD1.1 strAD1.1.1 ENST00000469289 ENSG00000243485
strAD1.1 strAD1.1.1 ENST00000469289 ENSG00000243485
答案2
也许只是:
< file cut -sd '"' -f2,4,8,10 | tr '"' ' '
即将输入视为"
分隔列的列表并提取第 2列、第 4列、第 8列和第 10列。
使用 GNU cut
,您可以将 替换| tr '"' ' '
为--output-delimiter=' '
.
这假设"
字符不会出现在行中的其他位置,那些gene_id
, transcript_id
... 属性始终出现并且始终按该顺序出现。
正如 Kusalananda 所指出的,您的示例中的情况并非如此,2,4,6,8
第一行和2,4,8,10
其他行应该是这样。
要进行更具表现力的匹配:仅应考虑第 9 个制表符分隔列并找到正确的属性名称,您可以使用正则表达式,例如:
< file pcregrep -o1 -o2 -o3 -o4 --om-separator=' ' '(?x)
^(?:[^\t]*+\t){8}(?=[^\t]*? \b gene_id \ +"([^"\t]*)")
(?=[^\t]*? \b transcript_id \ +"([^"\t]*)")
(?=[^\t]*? \b reference_id \ +"([^"\t]*)")
(?=[^\t]*? \b ref_gene_id \ +"([^"\t]*)")'
如果您没有pcregrep
或版本太旧而无法支持-o1...
,您可以使用perl
:
< file perl -lne 'print "$1 $2 $3 $4" if m{
^(?:[^\t]*+\t){8}(?=[^\t]*? \b gene_id \ +"([^"\t]*)")
(?=[^\t]*? \b transcript_id \ +"([^"\t]*)")
(?=[^\t]*? \b reference_id \ +"([^"\t]*)")
(?=[^\t]*? \b ref_gene_id \ +"([^"\t]*)")}x'
该正则表达式首先匹配前 8 个字段 ( (?:[^\t]*+\t){8}
),然后,我们有 4 个先行表达式 ( (?=...)
),因此我们将匹配这 8 个字段,前提是后面的内容匹配所有 4 个先行表达式。每个前瞻表达式都会查找属性之一并捕获值(在部件中(...)
)。这些捕获的值随后可在$1
、$2
、$3
、中使用$4
。
这样就允许属性以任何顺序排列。
请注意,它可能会被以下内容愚弄:
1 2 3 4 5 6 7 8 gene_id "transcript_id " ...
虽然可以解决这个问题,但可能不值得付出努力,因为我预计它不会出现在输入中。
当您使用 时,您还可以对第 9 个字段perl
进行更正式的解析。就像是:
< file perl -F'\t' -lane '
my %field;
while ($F[8] =~ /(\w+) +"(.*?)"/g) {$field{$1}=$2}
if (%field) {
print join " ", @field{
qw(gene_id transcript_id reference_id ref_gene_id
)}
}'
(这里,只要找到至少一个属性,就打印一行(与其他方法中请求的所有属性相反))。