从TCGA数据中提取lncRNA的表达量时,需要知道lncRNA的编号和对应的名称。这些信息可以从GTF文件中提取。
提取的话,可以采用如下的代码实现。
#!/usr/bin/perl -w
use strict;
my $biotype_file = shift @ARGV;
my $gtf = shift @ARGV;
my $biotype = shift @ARGV;
my %biotype_list;
open my $fh1, $biotype_file or die;
while (<$fh1>) {
chomp;
my @array = split /\t/, $_;
if($array[2]eq $biotype){
$biotype_list{$array[0]} = 1;
}
}
close $fh1;
open my $out, ">${biotype}_info.txt" or die;
print $out "Gene_id\tGene_id_info\tgene_name\tbiotype\n";
open my $fh2, $gtf or die;
while (<$fh2>) {
chomp;
next if /^#/;
my @array = split /\t/, $_;
next unless ($array[2] eq "gene");
$array[8] =~ /gene_id\s+"(\S+?)";.*gene_type\s+"(\S+?)";.*gene_name\s+"(\S+?)";/;
my $geneid = $1;
my $genebiotype = $2;
my $genename = $3;
my $gene_id_norm=(split("\\.",$geneid))[0];
if ($biotype_list{$genebiotype}) {
print $out "$gene_id_norm\t$geneid\t$genename\t$genebiotype\n";
}
}
close $fh2;
如果您对TCGA数据挖掘感兴趣,请学习我们的TCGA相关课程:
《GSEA富集分析》
如果觉得我的文章对您有用,请随意打赏。你的支持将鼓励我继续创作!