2022-04-21 用VDJdb数据库注释抗原表位、VDJmatch

目标图如下:

>>>开始

先从用VDJdb开始,上传了vdjtools格式、未经mixcr转换格式的文件,表头是:


在VDJdb官网不知为何注释失败,在网上找VDJdb使用教程的时候找到了用immunarch包进行注释的教程(使用immunarch包进行单细胞免疫组库数据分析(九):Annotate clonotypes using immune receptor databases - 简书 (jianshu.com)
)。


immunarch官网上用VDJdb数据库进行注释的步骤

因为groovy安装有点问题,安装了windows的installer,配置了环境变量,但是在命令行里groovy -v是没有的,故不能根据上图步骤执行。
下载的“vdjdb-db-master”数据库文件夹里没有代码所需的注释文件(vdjdb.slim.txt),到github上找,后来根据官网上的链接(”https://gitlab.com/immunomind/immunarch/-/blob/master/private)找到并下载,但是缺失了最后两列“clain”、“Pathology”,所以不能筛选疾病类型。

vdjdb =read.table("D:\\vdjdb.slim.txt",sep='\t',header=T,quote = "", fill=TRUE)
#, sep='\t',header=F
fix(vdjdb)
dim(vdjdb)
#colnames(vdjdb)=vdjdb[1,]
#rownames(vdjdb)=vdjdb[,1]
#vdjdb<-vdjdb[-1,]
colnames(data$A1.TRA)
colnames(immdata$data)
rownames(vdjdb)
vdjdb
vdjdb_st =read_tsv("D:\\SearchTable-2019-10-17 12_36_11.989.tsv","vdjdb-search", .species = "HomoSapiens", .chain = "TRA",.pathology = "ccRCC")
fix(vdjdb_st)
vdjdb_st
colnames(data)
dim(data)
colnames(vdjdb_st)=vdjdb_st[1,]
vdjdb_st<-vdjdb_st[-1,]
View(data)
library(immunarch)
data(immdata)
Convert.A1.TRA<-repLoad("C:\\Users\\Administrator.DESKTOP-4UQ3Q0K\\Downloads\\vdjtools-1.2.1\\vdjtools-1.2.1\\弦图Convert\\ConvertA1TRA.txt")
Convert.A1.TRA<-repLoad("D:\\mixcr-3.0.13\\弦图\\A1.TRA.txt")

Convert.A2.TRA<-repLoad("C:\\Users\\Administrator.DESKTOP-4UQ3Q0K\\Downloads\\vdjtools-1.2.1\\vdjtools-1.2.1\\弦图Convert\\Convert.A2.TRA.txt")

file_path = "D:\\mixcr-3.0.13\\弦图\\N.A"
A1.TRA<- repLoad(file_path)
View(A1.TRA)
A1.TRA

#Convert.A1.TRA<-as.list(Convert.A1.TRA)
#Convert.A2.TRA<-as.list(Convert.A2.TRA)
#sapply(data, function(element) {print(names(element))})
#data<-list(Convert.A1.TRA,Convert.A2.TRA)
names(data)[1] <- "A1.TRA"
names(data)[2] <- "A2.TRA"
class(data)
dbAnnotate(A1.TRA$data,vdjdb,"CDR3.aa", "cdr3")
dbAnnotate(immdata$data, vdjdb, "CDR3.aa", "cdr3")
dbAnnotate(immdata$data, vdjdb_st, "CDR3.aa", "CDR3")

最后一行的注释代码里,"cdr3aa"是自己的数据里的。"cdr3"是数据库的。
更改读取数据方式为repLoad以后,读入的数据类型就是list,不需要转换。(之前用read.table的方式读入,需要很麻烦的转换成二层嵌套的list结构,但其实不应该这样)
读入不了数据(横线处:unsupported format,skipping),应该是不能用转换格式后的数据。

改成读mixcr直接出来的数据,放在“D:\mixcr-3.0.13\弦图\N.A”(A1.TRA-A6.TRA)。还需要上传一个metadata.txt,里面要有一列列名为Sample,每行为样品名的矩阵,样品名可直接用数字简单代替。
得到用immunarch注释的结果如下:

但这并不是一开始想做的表格。重新思考发现,应该是用VDJtools做的表格,其中注释的命令下这个文件的后缀比较像,但是里面的列名内容不太一样。

官网写已弃用,抱着先试试的心态用命令行跑一下:

VDJtools命令行告知到vdjdb.cdr3.net这个网站,但这个网站就是一开始注释不成功的网站。接下来尝试用VDJmatch,先找到了github上的内容,检索发现也许是能产生预期的表格的。



>>>开始VDJmatch

VDJmatch网站:antigenomics/vdjmatch: ⚙️ Matching T-cell repertoire against a database of TCR antigen specificities (github.com)

1.3.1的版本不能用,可以下载旧一点的版本。下载了1.2.2的版本。
下载成功后,切换到所在文件夹,更新:java -jar vdjmatch-1.2.2.jar Update(不知道为什么不行,但是好像不影响之后的注释)

仿造运行vdjtools的命令输入java -jar vdjmatch-1.2.2.jar VDJmatch V1.2.2

试着用match一下看看能不能注释:java -jar ./vdjmatch-1.2.2.jar match A1.TRA.txt A1.TRA

报Missing required options: S, R
一开始不知道是什么意思,在github上VDJmatch 命令行选项往下看到:

根据VDJtools上填写的方式,参数可能填在功能名后面、样品名前面:java -jar ./vdjmatch-1.2.2.jar match -S human -R TRA A1.TRA.txt A1.TRA

看起来命令可以,数据库文件找不到,把之前下载的vdjdb.slim.txt复制到目录下,并修改文件名为vdjdb.slim.meta.txt

报找不到vdjdb.slim.txt
看样子程序需要两个文本文件
本来想和找vdjdb.slim.txt一样如法炮制到github上找vdjdb.slim.meta.txt,突然发现vdjmatch-1.2.2文件夹下有一个压缩的vdjdb文件夹,解压之后,产生了:

于是vdjdb.slim.meta.txt就这样送到了面前。(好奇驱使我先打开html看看是什么,发现打不开)
重复运行刚才的命令:

出来了一个A1.TRA.annot.summary文档,有sample_id metadata_blank db.column.name db.column.value unique frequency reads db.unique weight 但是是空的。

输入文件是这样的:

我认为也许用转换格式后的文件java -jar ./vdjmatch-1.2.2.jar match -S human -R TRA Convert.A1.TRA.txt A1.TRA
奇迹发生啦

Problem solving!

产生的两个文件

A1.TRB注释遇到内存不够的情况

但是A2.TRA很快注释完了,查看原文件发现文件大小是递增的,所以与原文件大小无关。

A3.TRA也是很快注释完了,先把TRA的都运行了一遍,都很快注释完成。
也许是TRB与TRA不同的原因,运行A2.TRB和A3.TRB时,都遇到java heap space,是时候再看看它还有哪些原因了。
先尝试了一下上次解决同样问题时用到的export MAX_MEMORY_OVERRIDE=12800
但是这个命令是liunx的,下面寻找在windows上设置的办法。

解决TRB的问题:

linux上VDJmatch的尝试

因为linux上运行快,并且也许不会遇到内存限制的问题。将文件复制到目录下。

java -jar vdjmatch-1.2.2.jar match -S human -R TRB Convert.A1.TRB.txt A1.TRB

报: java.net.ConnectException: Connection refused (Connection refused)
试了下TRA也是同样报错,所以应该不是A和B的关系。

先不用linux,换个角度,解决windows上内存不够的情况。在网上找到的方法:


>>>在windows上设置增加运行内存

新建环境变量行不通,搜索MAVEN _ OPTS

搜索maven环境搭建

要下载对应的包,虽然有点疑虑,因为TRA不用这样,但还是决定先试试。这个教程里写的很清楚。
Maven 环境配置 | 菜鸟教程 (runoob.com)
在这里还看到了mvn -v是linux和mac的命令。
两种方法(下图)。但依然报'mvn' 不是内部或外部命令,也不是可运行的程序

问题找到了,配置path变量的时候前面多打了一个分号,并且要重新打开命令行进行测试是否配置成功mvn -v或者mvn -version都可以。
重新跑TRB的命令,发现切换路径以后,直接
java -jar ./vdjmatch-1.2.2.jar match -S human -R TRB Convert.A1.TRB.txt A1.TRB就可以,不用像VDJtools一样先加载。
依然是Java heap space。
set MAVEN_OPTS=”-Xmx20000M -XX:MaxPermSize=512m”

依旧不行。根据上次的一次尝试,将用户变量里的MAVEN_OPTS如下更改:
-Xms4000m -Xmx20000m -XX:MaxPermSize=1024m -XX:+CMSClassUnloadingEnabled XX:SurvivorRatio=8 -XX:+UseConcMarkSweepGC -XX:+UseParNewGC -XX:+UseCMSCompactAtFullCollection -XX:+CMSParallelRemarkEnabled -XX:CMSInitiatingOccupancyFraction=70

老师之前成功的设置是在linux上用export MAX_MEMORY_OVERRIDE=12800,根据网络上其中一种在linux上的解决方法是export MAVEN_OPTS=-Xmx1024m -XX:MaxPermSize=512m,已知windows上设置用户变量的方式是把MAVEN_OPTS设置成-Xms4000m -Xmx20000m -XX:MaxPermSize=1024m,所以尝试在windows里把MAX_MEMORY_OVERRIDE设置为12800(猜测这个办法应该不行,在搜MAX_MEMORY_OVERRIDE的时候找到了另一种方法windows下配置tomcat服务器的jvm内存大小的两种方式 - 海绵般汲取 - 博客园 (cnblogs.com))明天试一下用tomcatw
再尝试:set "JAVA_OPTS=%JAVA_OPTS% -Xms4000m -Xmx20000M"(不行)
java -Xms4000m -Xmx20000 vdjmatch

查看堆的初始值和最大值
java -XX:+PrintFlagsFinal -version | findstr /i "HeapSize PermSize ThreadStackSize"

windows下查看某个java进程运行所需内存

C:\Program Files (x86)\Java\jdk1.8.0_73\bin下找到了jconsole

窗口看不到PID号,怎么都缩放不了。
bin下找到了jps,但是打开不了,会闪退。jstat也是。
更改了电脑的高级缩放设置为200%(之前为300%),终于看到了PID号

选中vdjmatch并点连接

已知在cmd控制台下输入 jps 命令,即可列出当前电脑运行的java程序的所有进程,但是

有写入权限,依然无法用jps,教程最后都是打开这个文件夹的图,不知道是不是正在运行的java进程的pid号的意思。(如果是的话我的如下图是7824)

教程里一般是发现“组或用户名中没有我当前使用的用户”,但是我的组或用户名中是有我的账户的,并且也可以在这个文件夹下创建文件。


>>>“后数据处理”

师姐想把前六个为一组进行注释,所以我用R将它们按列拼接了起来,拼接后还是vdjtools格式,但是程序就识别不了了(从count开始识别不了)。老师说多个样品的数据应该不能那么简单的拼接,应将6个注释结果汇总。

我发现拼接完后缺失了一些数据,所以接下来按注释后frequency这一列的值进行排序,把最优的筛选出来。


遇到的一些报错解决:

传入txt文件时读取不了的报错Error in type.convert.default data ill, as.is s= as.is[i],dec = dec, invalid multibyte string at '<ff><fe><63>'


错误保存方式示范

正确保存方式示范
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 214,922评论 6 497
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 91,591评论 3 389
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 160,546评论 0 350
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 57,467评论 1 288
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 66,553评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,580评论 1 293
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,588评论 3 414
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,334评论 0 270
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,780评论 1 307
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,092评论 2 330
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,270评论 1 344
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,925评论 5 338
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,573评论 3 322
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,194评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,437评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,154评论 2 366
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,127评论 2 352

推荐阅读更多精彩内容