007_10X突变位点分析和infercnv对比（一）

最新推荐文章于 2024-10-10 07:26:18 发布

TTS56

最新推荐文章于 2024-10-10 07:26:18 发布

阅读量915

点赞数

分类专栏：单细胞

本文链接：https://blog.csdn.net/qq_40966210/article/details/114016140

版权

单细胞专栏收录该内容

19 篇文章 21 订阅

订阅专栏

VarTrix 突变位点分析对比inferCNV（一）

动机：记录10x公司开发的软件工具VarTrix和inferCNV，在肿瘤细胞识别中的区别和联系

我一开始以为这个很简单，一天就能弄完，没想到后来问题越来越多，时间越拖越久。。

我觉得我离被拉入黑名单不远了

Cellranger 流程

下载10X的测试数据

因为主要是针对肿瘤细胞的，于是下载10x公司的胶质瘤测试数据

# 下载10x的胶质瘤数据，主要是为了测试一下infercnv和这个10x开发的软件的具体区别在哪里
wget -c https://cf.10xgenomics.com/samples/cell-exp/4.0.0/Parent_SC3v3_Human_Glioblastoma/Parent_SC3v3_Human_Glioblastoma_fastqs.tar

mv /mnt/d/迅雷下载/Parent_SC3v3_Human_Glioblastoma_fastqs.tar /home/yuansh/10XtestData/
# 看一下下载的数据是不是和官网的一样
md5sum Parent_SC3v3_Human_Glioblastoma_fastqs.tar
# 解压（我使用的是zsh编译，如果是bash的话要使用tar -vzxf Parent_SC3v3_Human_Glioblastoma_fastqs.tar）
x Parent_SC3v3_Human_Glioblastoma_fastqs.tar 
# 解压好了就删掉就行
# rm *.tar

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-pGm83Dhp-1614139374962)(…/…/…/image/image-20201211162431120.png)]

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-gOsMZOim-1614139374965)(…/…/…/image/image-20201211162442359.png)]

然后跑一下cellranger流程

注意事项：

ref是参考基因组，10x官网下载的
fastqFilePath 这个是解压完的文件夹的路径
id是 _S1前面的东西，不包括
localcores 是指电脑核心个数，如果是6和就设置6就行了，不要设置太多会报错，但是好像有的电脑设置超出了也不会报错，反正就按照实际情况设置，如果你的电脑最高核心数是6，那么设置4是一定不会报错的。然后如果要把核心占满尽量把后台所有的软件全部关闭。否则即便设置满核心照样没什么用处

mkdir cellranger_result && cd cellranger_result
ref=/home/yuansh/bioreference/10x/grch38_10x
fastqFilePath=/home/yuansh/10XtestData/Parent_SC3v3_Human_Glioblastoma_fastqs/Parent_SC3v3_Human_Glioblastoma
id=Parent_SC3v3_Human_Glioblastoma 
# 通常情况下应该是大于1个小时的
cellranger count --id=$id \
                  --transcriptome=$ref \
                  --fastqs=$fastqFilePath \
                  --sample=$id \
                  --localcores=15\
                  --nosecondary \
                  --localmem=30

VarTrix流程

首先看一下官方对 VarTrix 软件的介绍

VarTrix is a software tool for extracting single cell variant information from 10x Genomics single cell data. VarTrix will take a set of previously defined variant calls and use that to identify those variants in the single cell data. VarTrix does not perform variant calling. VarTrix is useful for evaluating heterogeneity within a sample, which means that the types of variants that will be useful are either somatic or contained within a copy number variant (CNV) event. which means that the types of variants that will be useful are either somatic or contained within a copy number variant (CNV) event .

意思就是，他这个结果，不仅是突变位点的信息，还包括拷贝数变异的信息，这就是为什么我要把它和infercnv进行比较

软件下载

https://github.com/10XGenomics/vartrix/releases

下载最新版本的就行

在这里插入图片描述

然后直接解压就可以直接用

使用方法

最主要的就是vcf文件的获取，刚开始没有注意到导致跑了三四遍结果都是奇奇怪怪的

注意到官方文档的一句话:

Pre-called variants to be used as input to VarTrix can be generated in many different ways such as gathering calls from existing variant databases or performing variant calling on bulk or single cell genome or transcriptome data. It is important to note that generating variants from bulk or single cell RNA-seq datasets is challenging. Noise inherent in reverse transcription leads to a high false positive rate . We recommend looking at the Broad Institute’s GATK and Mutect2 best practices guide for calling variants in RNAseq. An alternative approach is to determine somatic variants using WGS data generated from the same sample as the scRNA-seq library.

~~意思应该是走一遍GATK流程生成的vcf才能用~~

然而并不是这个意思，10x官网有给出一个全基因组的calling代码，我这个是转录组的数据，然后我强行走了一遍，没有报错，但是最后的结果太奇怪了。（因为没有保存，所以无法展示）

讲实话我其实完全没看懂上面那一段话到底什么意思。（我觉得并不是我英语水平的问题）

上文中提到，可以使用基因组或者转录组进行calling，但是I0x官方的SN calling代码没有转录组的代码。然后按照上面的意思应该是随便找一个差不多的bulk数据的vcf文件也行，然后我下载了gatk的谷歌云hg38的confidence.vcf文件和omin.vcf文件进行vartrix结果都是错误的，于是我干脆把所有的hg38涉及到的vcf全部下载下来。写了一个for循环直接看到底那一个可以用，最后的结果就是：

dbsnp_146.hg38.vcf.gz 这个文件是可以用的

但还是我还是很想知道performing variant calling on bulk or single cell genome or transcriptome data.这个到底是什么意思，如果转录组可以的话，要如何calling

#vartrix \ 你们调用的时候注意一下，我是因为把名字给改了
#-v <path_to_input_vcf> \ 自己准备的vcf文件
#-b <path_to_cellranger_bam> \ cellranger输出的单细胞bam文件
#-f <path_to_fasta_file> \ 自己准备的vcf文件
#-c <path_to_cell_barcodes_file> \ cellranger 输出的文件
#--threads 30 \ 设置线程
#-o <path_for_output_matrix> 输出文件路径


# 这个流程需要太多的时间了，至少12个小时（我这还仅仅只是一个样本的数据）
# vcf=/home/yuansh/bioreference/gatk/dbsnp_146.hg38.vcf.gz
# bam=/home/yuansh/10XtestData/Parent_SC3v3_Human_Glioblastoma_fastqs/cellranger_result/Parent_SC3v3_Human_Glioblastoma/outs/possorted_genome_bam.bam
# fa=/home/yuansh/bioreference/10x/grch38_10x/fasta/genome.fa
# tsv=/home/yuansh/10XtestData/Parent_SC3v3_Human_Glioblastoma_fastqs/cellranger_result/Parent_SC3v3_Human_Glioblastoma/outs/filtered_feature_bc_matrix/barcodes.tsv.gz
# 
# vartrix \
# -v $vcf \
# -b $bam \
# -f $fa \
# -c $tsv \
# --threads 32 \
# -o out/VarTrix.mtx
	
#根据上文官方文档的意思，应该是可以直接使用单细胞的数据进行call出vcf文件的（虽然我之前跑过一次，并且结果是不对的，我这次打算在跑一次）
bcftools mpileup -Ou -f /home/yuansh/bioreference/10x/grch38_10x/fasta/genome.fa \
/home/yuansh/10XtestData/Parent_SC3v3_Human_Glioblastoma_fastqs/cellranger_result/Parent_SC3v3_Human_Glioblastoma/outs/possorted_genome_bam.bam | bcftools call -Ou -mv | bcftools filter -s LowQual -e '%QUAL<20 || DP>100' > my.vcf

vcf=/home/yuansh/my.vcf
bam=/home/yuansh/10XtestData/Parent_SC3v3_Human_Glioblastoma_fastqs/cellranger_result/Parent_SC3v3_Human_Glioblastoma/outs/possorted_genome_bam.bam
fa=/home/yuansh/bioreference/10x/grch38_10x/fasta/genome.fa
tsv=/home/yuansh/10XtestData/Parent_SC3v3_Human_Glioblastoma_fastqs/cellranger_result/Parent_SC3v3_Human_Glioblastoma/outs/filtered_feature_bc_matrix/barcodes.tsv.gz

vartrix \
-v $vcf \
-b $bam \
-f $fa \
-c $tsv \
--threads 32 \
-o out/VarTrix.mtx

# 将vcf转换为seurat可读文件
# 转换为seurat可读对象
awk '{print $1,$2}' my.vcf > out/SNV.loci.txt
sed -i 's/\s/:/g' out/SNV.loci.txt 
head -n 50 |tail -n30 out/SNV.loci.txt
# 将文件移动至桌面
cp -r out/ /mnt/c/Users/yuansh/Desktop/VartrixVsInferCNV

要点：

vcf文件必须和上面的cellranger是一样的版本，比如我上面cellranger流程使用的是grch38，则这里就要用grch38对应的vcf文件
vcf文件选择：根据官方给的单细胞vcf应该可以作为输入文件但是不知道如何操作，dbsnp_146.hg38.vcf 这个文件是一定能用的
hg38的fa文件下面必须要有配套的索引文件 .fai文件
bam文件用的是filter中的bam

基本上调用了全部的核，1个样本大概需要10个小时（主要还是和这个参考vcf有关，如果使用单细胞call出来的vcf只需要4个小时）

下面必须要有配套的索引文件 .fai文件