分类目录归档:Default Category

统计GTF文件中转录本的长度 Calculate transcript length from gtf file

gtf 文件INPUT

chr1    PacBio  exon    763020  763155  .       +       .       gene_id "XLOC_000001"; transcript_id "TCONS_00000001"; exon_number "1"; gene_name "LINC01128"; oId "PB.5.1"; nearest_ref "NR_047519"; class_code "j"; tss_id "TSS1";
chr1    PacBio  exon    764383  764484  .       +       .       gene_id "XLOC_000001"; transcript_id "TCONS_00000001"; exon_number "2"; gene_name "LINC01128"; oId "PB.5.1"; nearest_ref "NR_047519"; class_code "j"; tss_id "TSS1";
chr1    PacBio  exon    776580  776753  .       +       .       gene_id "XLOC_000001"; transcript_id "TCONS_00000001"; exon_number "3"; gene_name "LINC01128"; oId "PB.5.1"; nearest_ref "NR_047519"; class_code "j"; tss_id "TSS1";

Code

awk -F"\t" '
{
if ($3=="exon")
    {
        split($9, a, " ");
        ID=a[4]; 
        L[ID]+=$5-$4+1
    } 
}
END{
    for(i in L)
        {print i"\t"L[i]}
    }
' gtf-file >output

You can change array index and feature type to choose what you want to count. For example, if you want to calculate gene length, change exon to gene (make sure your gtf file has this feature).

输出OUTPUT

"TCONS_00052437";       5950
"TCONS_00049398";       988
"TCONS_00031225";       6005
"TCONS_00026369";       825

参考:https://gist.github.com/sp00nman/10372555

RefSeq的gtf文件

注释有很多版本,比如ensembl,gencode, ucsc known gene, NCBI的RefSeqGene。最近就需要NM id的注释,但NCBI提供的是gff3格式的,而且很乱。用UCSC table browser下载的gtf版本的RefSeq,没有转录本和基因之间的关系,也没有基因symbol。

比如Ensembl,其实Ensembl的gtf挺好用的,不过这次我因为需要NM编号的注释(笨方法是将ensembl id转成NCBI的refSeq的ID,但这不是最优的方法,ID mapping有可能对不上,不如直接用NM的注释)。

chr1    ensembl_havana  gene    11869   14412   .       +       .       gene_id "ENSG00000223972"; gene_version "4"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene";
chr1    havana  transcript      11869   14409   .       +       .       gene_id "ENSG00000223972"; gene_version "4"; transcript_id "ENST00000456328"; transcript_version "2"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; transcript_name "DDX11L1-002"; transcript_source "havana"; transcript_biotype "processed_transcript"; havana_transcript "OTTHUMT00000362751"; havana_transcript_version "1"; tag "basic";
chr1    havana  exon    11869   12227   .       +       .       gene_id "ENSG00000223972"; gene_version "4"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "1"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; transcript_name "DDX11L1-002"; transcript_source "havana"; transcript_biotype "processed_transcript"; havana_transcript "OTTHUMT00000362751"; havana_transcript_version "1"; exon_id "ENSE00002234944"; exon_version "1"; tag "basic";

UCSC table browser下载的refGene的gtf,这个文件不对的地方是gene id和transcript id是一个,而我需要gene和transcript的关系

chr1    hg19_refGene    exon    66999276        66999355        0.000000        +       .       gene_id "NM_001308203"; transcript_id "NM_001308203";
chr1    hg19_refGene    start_codon     67000042        67000044        0.000000        +       .       gene_id "NM_001308203"; transcript_id "NM_001308203";
chr1    hg19_refGene    CDS     67000042        67000051        0.000000        +       0       gene_id "NM_001308203"; transcript_id "NM_001308203";
chr1    hg19_refGene    exon    66999929        67000051        0.000000        +       .       gene_id "NM_001308203"; transcript_id "NM_001308203";
chr1    hg19_refGene    CDS     67091530        67091593        0.000000        +       2       gene_id "NM_001308203"; transcript_id "NM_001308203";
chr1    hg19_refGene    exon    67091530        67091593        0.000000        +       .       gene_id "NM_001308203"; transcript_id "NM_001308203";

中间找了很多,终于找到一个方法,https://bioinformatics.stackexchange.com/questions/2548/hg38-gtf-file-with-refseq-annotations

mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -N -e "select * from refGene" hg19 | cut -f2- | genePredToGtf -source=hg19.refGene.ucsc file stdin stdout > refSeq.hg19.gtf

genePredToGtf可以从UCSC下载,需要注意的直接下载refFlat的文件(不通过mysql),用genePhredToGtf生成的文件和UCSC table browser下载的gtf文件是一样的,通过mysql下载用genephreadToGtf是我想要的。

比如这样,这才是我想要的

chrX    hg19.refGene.ucsc       exon    70683674        70685855        .       +       .       gene_id "TAF1"; transcript_id "NM_138923"; exon_number "38"; exon_id "NM_138923.38"; gene_name "TAF1";
chrX    hg19.refGene.ucsc       CDS     70683674        70683893        .       +       1       gene_id "TAF1"; transcript_id "NM_138923"; exon_number "38"; exon_id "NM_138923.38"; gene_name "TAF1";
chrX    hg19.refGene.ucsc       start_codon     70586225        70586227        .       +       0       gene_id "TAF1"; transcript_id "NM_138923"; exon_number "1"; exon_id "NM_138923.1"; gene_name "TAF1";
chrX    hg19.refGene.ucsc       stop_codon      70683894        70683896        .       +       0       gene_id "TAF1"; transcript_id "NM_138923"; exon_number "38"; exon_id "NM_138923.38"; gene_name "TAF1";

Ref: https://bioinformatics.stackexchange.com/questions/2548/hg38-gtf-file-with-refseq-annotations
#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################

TCGABiolinks下载TCGA数据做生存分析

以前的工作是全基因组或全外分析,不涉及癌症和生存分析,但现在的工作主要围绕癌症方面,生存分析一定少不了。实验室小伙伴推荐用TCGAbiolinks下载TCGA的数据,于是研究了如何用TCGABiolinks下载TCGA的数据,以下载RNA的count数据为例,并做生存分析。

# 安装相关的包
if (!requireNamespace("BiocManager", quietly=TRUE))   
  install.packages("BiocManager")
BiocManager::install("TCGAbiolinks")
BiocManager::install("SummarizedExperiment")
install.packages('survival')
install.packages('survminer')
library(SummarizedExperiment)
library(TCGAbiolinks)
library(survival)
library(survminer)

# 项目的概括信息,project的名称可以从 https://portal.gdc.cancer.gov/ 上面选择对应的器官,进入之后左侧列表中就会显示。我们以下载TCGA-LIHC项目的数据为例。
TCGAbiolinks:::getProjectSummary("TCGA-LIHC") 
# 获取该项目样本的临床信息
clinical <- GDCquery_clinic(project = "TCGA-LIHC", type = "clinical")

# 下载TCGA-LIHC项目的rna-seq的counts数据,构建GDCquery,具体的数据类型的写法可以参考 https://portal.gdc.cancer.gov/repository 左侧列表
# 比如data type有Raw Simple Somatic Mutation, Annotated Somatic Mutation, Aligned Reads, Gene Expression Quantification, Slide Image这几种,不过不同的项目,有可能有不同的数据类型。
query <- GDCquery(project = "TCGA-LIHC", 
		    experimental.strategy = "RNA-Seq",
                    data.category = "Transcriptome Profiling", 
                    data.type = "Gene Expression Quantification",
		    workflow.type = "HTSeq - Counts")
# 下载数据
GDCdownload(query)
# 导入数据
LIHCRnaseq <- GDCprepare(query)
# 得到样本的基因counts矩阵,每行为一个基因,每列为一个样本
count_matrix <- assay(LIHCRnaseq)

# TCGA样本的编号以'-'分割,前三列是患者编号,第四列是类型,一般11表示正常样本,01表示肿瘤样本,见https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables/sample-type-codes
samplesNT <- TCGAquery_SampleTypes(barcode = colnames(count_matrix),typesample = c("NT"))
# selection of tumor samples "TP"
samplesTP <- TCGAquery_SampleTypes(barcode = colnames(count_matrix),typesample = c("TP"))

# 根据ensembl id选择特定基因在癌症样本中的表达
gexp <- count_matrix[c("ENSG00000198431"),samplesTP]

# 这样选出的样本都是肿瘤样本,但样本名称是TCGA-CC-A7IL-01A-11R-A33R-07,与clinical信息不对应,需要将名字简化成TCGA-CC-A7IL
names(gexp) <-  sapply(strsplit(names(gexp),'-'),function(x) paste0(x[1:3],collapse="-"))

# 将表达数据和clinical数据合并,并且挑出来用于做生存分析的数据
clinical$"ENSG00000198431" <- gexp[clinical$submitter_id]
df<-subset(clinical,select =c(submitter_id,vital_status,days_to_death,ENSG00000198431))

#去掉NA
df <- df[!is.na(df$ENSG00000198431),]

#取基因的平均值,大于平均值的为H,小于为L
df$exp <- ''
df[df$ENSG00000198431 >= mean(df$ENSG00000198431),]$exp <- "H"
df[df$ENSG00000198431 <  mean(df$ENSG00000198431),]$exp <- "L"

# 将status表示患者结局,1表示删失,2表示死亡
df[df$vital_status=='Dead',]$vital_status <- 2
df[df$vital_status=='Alive',]$vital_status <- 1
df$vital_status <- as.numeric(df$vital_status)

# 建模
fit <- survfit(Surv(days_to_death, vital_status)~exp, data=df) # 根据表达建模
# 显示P value
surv_pvalue(fit)$pval.txt 
# 画图
ggsurvplot(fit,pval=TRUE)

继续阅读

误删hyper-v的avhdx文件

因为对hyper-v不是很熟悉,点了一下检查点,生成了一个avhdx文件,这个文件其实后续hyper-v会将其合并到vhdx的虚拟磁盘中。而我当时手贱手工的删除了avhdx文件,导致hyper-v找不到这个文件,vhdx也挂起等待合并,虚拟机迟迟不能启动。

有一种解决办法是文件恢复,但我用了几个文件都没有恢复成。实验室师兄(超级牛)新建了一个虚拟机挂载已有的vhdx文件,尝试用vhdx文件启动,显示不能启动,但在新的虚拟机下没有提示要合并,提示老系统的vhdx还有戏。

于是又新建了一个虚拟机实例,创建虚拟机实例之后,尝试将以前的vhdx文件挂载到新的虚拟机上,重启发现竟然以老的系统启动了。感谢能够启动,避免实验室的数据丢失。

根据结果反推,第一个shimx64.efi和Ubuntu.vhdx都是以前的系统,第二个shimx64.efi新的虚拟机的,硬盘驱动器已经换成了老系统的。

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################

Fusion Gene Annotation

STAR-FUSION和FusonAnnotator都属于Trinity Trinity Cancer Transcriptome Analysis Toolkit Fusion-finding modules。
CTAT_HumanFusionLib现阶段整合了各种资源帮助分析癌症生物学相关的fusion,同样也鉴别可能在正常样本只能出现的fusion。下载地址:https://data.broadinstitute.org/Trinity/CTAT_RESOURCE_LIB/

FusionAnnotator –genome_lib_dir GRCh37_gencode_v19_CTAT_lib_July192017/ctat_genome_lib_build_dir/ \
–annotate fusions.list.txt
fusions.list.txt为star-fusion的结果中的第一列,两个参与融合的基因中间用–连在一起,就可以用FusionAnnotator进行注释,相关的标签会注释到融合基因上。

会有三类标签,每类下面又有很多具体的来源标签:
Fusions relevant to cancer biology
Individual genes of cancer relevance, which may show up in fusions
Red Herrings: Fusion pairs that may not be relevant to cancer, and potential false positives.

通过注释,就可以了解到分析结果中的融合基因是否在其他数据库中出现过,或者可能是和癌症无关的突变。

参考:https://github.com/FusionAnnotator/CTAT_HumanFusionLib/wiki