目标图如下:
>>>开始
先从用VDJdb开始,上传了vdjtools格式、未经mixcr转换格式的文件,表头是:
在VDJdb官网不知为何注释失败,在网上找VDJdb使用教程的时候找到了用immunarch包进行注释的教程(使用immunarch包进行单细胞免疫组库数据分析(九):Annotate clonotypes using immune receptor databases - 简书 (jianshu.com)
)。
因为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
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>'