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

NIS+NFS+SGE

需求,把多台服务器组成一个cluster(SGE),把一台电脑(比如存储)的home文件件共享给其他服务器(NFS),共用一个home文件夹,并进行用户的统一管理(NIS)。

操作系统为操作系统:CentOS,用virtual box虚拟出来的系统做测试。
server端:10.0.2.5
client或compute端:在同样网段

1,NFS共享存储

通过nfs,实现每台服务器都有同样的路径和文件,便于后续集群管理。这里共享两个路径,一个是server端的/home路径,实现每个服务器都有同样的家目录,一个是/opt/gridengine用于安装SGE。

1.1 Server端:

安装相关软件,NFS的端口是不固定的(因此如果客户端连不上的时候,往往需要iptables -F清理一下),客户端要准确的获得NFS服务器所使用的端口,就需要RPC服务。RPC最主要的功能就是记录每个NFS功能所对应的端口号,并且在NFS客户端请求时将该端口和功能对应的信息传递给请求数据的NFS客户端,让客户端可以链接到正确的端口上去,从而实现数据传输。

yum install nfs-utils rpcbind

开机启动rpcbind

systemctl enable rpcbind.service

开启rpcbind

systemctl start rpcbind.service

设置要共享的目录

mkdir /opt/gridengine
vi /etc/exports
/home 10.0.2.0/255.255.255.0(rw,sync)
/opt/gridengine 10.0.2.0/255.255.255.0(rw,sync)

nfs开机启动和开启服务

systemctl enable nfs
systemctl start nfs

生效export

exportfs -r -v

1.2 客户端:

yum install nfs-utils rpcbind
systemctl enable  rpcbind.service
systemctl restart rpcbind.service
systemctl enable nfs
systemctl start  nfs

设置自动挂载,用tab分割,不是空格

mkdir /opt/gridengine
vi /etc/fstab
10.0.2.5:/home	/home	nfs	defaults	0	0
10.0.2.5:/opt/gridengine	/opt/gridengine	nfs	defaults	0	0

挂载

mount -a

这样的话,客户端服务器的home目录都是server端的home家目录。

2,NIS(Network Information Service)

通过NIS实现帐号的统一权限管理和认证,避免在多台服务器上重复开设帐号

2.1 Server端:

继续阅读