我有一个大的gtf
。我在这里分享eg.gtf
如下:
chr22 Cufflinks transcript 10695955 10696708 . + . transcript_id "first_11345700"; gene_id "XLOC_158970"; gene_name "XLOC_158970"; oId "TCONS_00353198"; class_code "u"; tss_id "TSS369767"; original_gene_id "XLOC_158970";
chr22 Cufflinks exon 10702915 10703826 . + . transcript_id "first_11345701"; gene_id "ENSG00000277248.1"; gene_name "ENSG00000277248.1"; exon_number "1"; original_gene_id "ENSG00000277248.1";
chr22 Cufflinks transcript 10702915 10707278 . + . transcript_id "first_11345701"; gene_id "ENSG00000277248.1"; gene_name "ENSG00000277248.1"; oId "TCONS_00353199"; class_code "u"; tss_id "TSS369769"; original_gene_id "ENSG00000277248.1";
我awk
使用列号提取了所需的字段,如下所示:
cat eg.gtf | awk 'OFS="\t" {if ($3=="transcript") {print $1,$4-1,$5,$12,$7}}' | tr -d '";'
输出如下所示:
chr22 10695954 10696708 XLOC_158970 +
chr22 10702914 10707278 ENSG00000277248.1 +
12
我不想在命令中使用列号awk
,而是想提取带有名称的字段。
笔记:12th column has different names, starting with E or X or M or N or S
。
12
如何在不给出命令中列号的情况下获取第 12 个字段awk
?有没有办法使用gene_id
第 11 列的术语获取 12 字段?
答案1
在每个 Unix 机器上的任何 shell 中使用任何 awk:
$ cat tst.awk
BEGIN {
FS=OFS="\t"
}
$3 == "transcript" {
n = split($NF,tmp,/[; "]+/)
for ( i=1; i<n; i+=2 ) {
vals[tmp[i]] = tmp[i+1]
}
print $1, $4-1, $5, vals["gene_id"], $7
}
$ awk -f tst.awk file
chr22 10695954 10696708 XLOC_158970 +
chr22 10702914 10707278 ENSG00000277248.1 +
答案2
GTF 是制表符分隔格式,您不能依赖子字段的顺序。您需要改用文本解析。我会用 perl 来做:
$ perl -F'\t' -lane '
if($F[2] eq "transcript"){
$gene_id = /gene_id\s*"([^"]+)"/ ? $1 : "None";
print join("\t",$F[0],$F[3]-1,$F[4],$gene_id,$F[6])
}' file
chr22 10695954 10696708 XLOC_158970 +
chr22 10702914 10707278 ENSG00000277248.1 +
这具有在没有基因名称的情况下打印“None”的优点(例如,非编码转录本可能没有基因名称,或者未知转录本等)。
答案3
使用磨坊主( ) 同时假设输入是制表符分隔的文件,其中第 9 个字段由- 分隔的键值字段mlr
组成(键用双引号值字符串中的空格分隔):;
mlr --inidx --ifs tab --otsv --headerless-csv-output \
filter '$3 == "transcript"' then \
put '$4 -= 1; $9 = gsub($9, "; ", ";"); $9 = gsub($9, "\"", "")' then \
nest --explode --pairs --across-fields --nested-fs ';' --nested-ps ' ' -f 9 then \
cut -o -f 1,4,5,gene_id,7 file.gtf
这会将文件读取为整数索引制表符分隔字段,并创建无标头制表符分隔数据输出 (TSV)。我们无法将数据读取为 TSV,因为 Miller 抱怨嵌入的引号(如果字段本身未加引号,则不能在字段中嵌入引号)。
处理从filter
丢弃第三个字段不是 string 的任何记录的操作开始transcript
。
然后,我们使用put
来调整第四个字段的值,将其减一。我们还删除了;
第 9 个字段(带有额外注释的字段)中每个字段后面的空格以及所有字段的双引号。
该nest --explode
操作创建新命名的通过“分解”其键值对,从第 9 个字段中的数据中提取字段。
然后我们可以使用cut
提取我们想要的字段,其中包括该gene_id
字段。
输出问题中的数据:
chr22 10695954 10696708 XLOC_158970 +
chr22 10702914 10707278 ENSG00000277248.1 +