融合基因鉴定以及GATK寻找突变

生信菜鸟团

发布于 2023-09-08 12:28:23

1.6K1

发布于 2023-09-08 12:28:23

文章被收录于专栏：生信菜鸟团

上周的癌症样本全转录组数据的融合基因鉴定中我们拿到数据进行一系列比对过滤后使用star完成了基因组比对，并通过设置参数拿到了Chimeric.out.junction文件以便star-fusion进行融合基因的鉴定

但是这后面遇到了一些问题，本周我们就着力解决，完成star-fusion鉴定融合基因，并使用gatk寻找突变

参考：RNA-seq数据分析完全指北-09：寻找融合基因

寻找融合基因

STAR-Fusion鉴定融合基因

STAR和STAR-Fusion的版本一定要相对应

参考：

融合检测starfusion的运用和版本问题（https://zhuanlan.zhihu.com/p/587614424）融合基因分析-STAR-Fusion（https://evvail.com/2022/10/05/2822.html） STAR Fusion release and CTAT Genome Lib Compatibility Matrix（https://github.com/STAR-Fusion/STAR-Fusion/wiki/STAR-Fusion-release-and-CTAT-Genome-Lib-Compatibility-Matrix）

使用 mamba search 检查star和star-fusion的版本，发现使用的star为最新版 2.7.10b，在表中无对应star-fusion版本，直接也下载star-fusion最新版 star-fusion=1.12.0（由于环境软件版本问题这里单独创建了一个python3.10的环境安装，在安装的时候发现自动安装了一个2.7.8a版本的star）

此外，还需安装相关数据库

参考：

融合基因分析-STAR-Fusion（https://evvail.com/2022/10/05/2822.html） https://data.broadinstitute.org/Trinity/CTAT_RESOURCE_LIB/ 其中’plug-n-play’是已经建好数据库的可以开箱直接使用，但是文件比较大下载需要时间;另一个是source版本的大约4G左右，下载好后需要自己建索引和必要的库文件融合基因分析工具：STAR-Fusion (https://zhuanlan.zhihu.com/p/586371542)

我们这里直接下载开箱即用的数据库，并解压https://data.broadinstitute.org/Trinity/CTAT_RESOURCE_LIB/GRCh38_gencode_v37_CTAT_lib_Mar012021.plug-n-play.tar.gz

cat srr_acc_list.txt | while read id
do 
        STAR-Fusion --genome_lib_dir ~/refdata/index/star_fusion/GRCh38_gencode_v37_CTAT_lib_Mar012021.plug-n-play/ctat_genome_lib_build_dir -J 4.mapping/${id}Chimeric.out.junction --output_dir 5.fusion/${id} --CPU 16
done

nohup bash fusion.sh &> fusion_log &

完成后查看：

$grep -w "complete" fusion_log 
 * STAR-Fusion complete.  See output: /home/data/t120455/fusion_and_mut/5.fusion/SRR11178348/star-fusion.fusion_predictions.tsv (or .abridged.tsv version)
 * STAR-Fusion complete.  See output: /home/data/t120455/fusion_and_mut/5.fusion/SRR11178349/star-fusion.fusion_predictions.tsv (or .abridged.tsv version)
 * STAR-Fusion complete.  See output: /home/data/t120455/fusion_and_mut/5.fusion/SRR11178350/star-fusion.fusion_predictions.tsv (or .abridged.tsv version)
 * STAR-Fusion complete.  See output: /home/data/t120455/fusion_and_mut/5.fusion/SRR11178351/star-fusion.fusion_predictions.tsv (or .abridged.tsv version)
 * STAR-Fusion complete.  See output: /home/data/t120455/fusion_and_mut/5.fusion/SRR11178352/star-fusion.fusion_predictions.tsv (or .abridged.tsv version)
 * STAR-Fusion complete.  See output: /home/data/t120455/fusion_and_mut/5.fusion/SRR11178353/star-fusion.fusion_predictions.tsv (or .abridged.tsv version)

未找到任何融合基因事件，只有表头信息

#FusionName     JunctionReadCount       SpanningFragCount       est_J   est_S   SpliceType      LeftGene
        LeftBreakpoint  RightGene       RightBreakpoint JunctionReads   SpanningFrags   LargeAnchorSupport
      FFPM    LeftBreakDinuc  LeftBreakEntropy        RightBreakDinuc RightBreakEntropy       annots

而原推文中是存在的，并且这是胃癌样本，按理来说也应该有基因融合事件

就这个问题我询问了曾老师

Q：胃癌的样本，原推文找出了融合基因事件，我跟着走下来没有找到，在想新版的软件、参考基因组和注释文件会有什么影响吗 A：如果是全新的参考基因组什么的，应该不会影响它。融合事件是癌症特有的，只要你是癌症样品，按照道理它就应该有。有没有可能是融合事件过滤的参数太严格了呢？

在前面star比对获得Chimeric.out.junction文件时，我们像原推文那样设置了一系列参数：

融合基因部分

--chimSegmentMin 20    # 每对嵌合读段较短的一端至少有20个碱基，即PE150允许280+20结构的融合基因
--chimOutJunctionFormat 1    # 输出的Chimeric.out.junction文件可直接用于融合基因
--chimSegmentReadGapMax 0    # 嵌合读段之间不允许空位
--chimJunctionOverhangMin 20    # 嵌合的junction的最少的overhang是20个碱基，为了过滤非常短的外显子，即连续剪切事件

于是我们这里更换star-fusion模式看看，直接对fq文件进行处理，跳过这些我们设置的参数

参考：

融合基因分析-STAR-Fusion（https://evvail.com/2022/10/05/2822.html）

融合基因分析工具：STAR-Fusion (https://zhuanlan.zhihu.com/p/586371542)

直接从双端按测序的fq文件开始：

cat srr_acc_list.txt | while read id
do 
        STAR-Fusion --genome_lib_dir ~/refdata/index/star_fusion/GRCh38_gencode_v37_CTAT_lib_Mar012021.plug-n-play/ctat_genome_lib_build_dir --left_fq 3.trimg/${id}_rmr_rmOther_1_val_1.fq.gz --right_fq 3.trimg/${id}_rmr_rmOther_2_val_2.fq.gz --output_dir 5.fusion_from_fq/${id} --CPU 16
done

完成后查看：

$grep -w "complete" fusion_from_fq_log 
 * STAR-Fusion complete.  See output: /home/data/t120455/fusion_and_mut/5.fusion_from_fq/SRR11178348/star-fusion.fusion_predictions.tsv (or .abridged.tsv version)
 * STAR-Fusion complete.  See output: /home/data/t120455/fusion_and_mut/5.fusion_from_fq/SRR11178349/star-fusion.fusion_predictions.tsv (or .abridged.tsv version)
 * STAR-Fusion complete.  See output: /home/data/t120455/fusion_and_mut/5.fusion_from_fq/SRR11178350/star-fusion.fusion_predictions.tsv (or .abridged.tsv version)
 * STAR-Fusion complete.  See output: /home/data/t120455/fusion_and_mut/5.fusion_from_fq/SRR11178351/star-fusion.fusion_predictions.tsv (or .abridged.tsv version)
 * STAR-Fusion complete.  See output: /home/data/t120455/fusion_and_mut/5.fusion_from_fq/SRR11178352/star-fusion.fusion_predictions.tsv (or .abridged.tsv version)
 * STAR-Fusion complete.  See output: /home/data/t120455/fusion_and_mut/5.fusion_from_fq/SRR11178353/star-fusion.fusion_predictions.tsv (or .abridged.tsv version)

直接从fq文件使用默认参数star-fusion找到了融合基因，以SRR11178349样本为例：

/fusion_and_mut/5.fusion_from_fq/SRR11178348/star-fusion.fusion_predictions.tsv

$grep -v "^#" /home/data/t120455/fusion_and_mut/5.fusion_from_fq/SRR11178348/star-fusion.fusion_predictions.tsv | awk '{print $1}'
E2F3--CDKAL1
WNK1--ERC1
PHACTR4--RCC1
GAS6--TUBB2B
AC022400.4--BMS1P4
WWTR1--COMMD2
TMEM181--WTAP
AL627309.5--SEPTIN14
MED13--DDX46
LINC02246--NRIP1
LSM14A--GARRE1
ERICH2--AC007405.3

解读结果文件：

参考：

融合基因分析工具：STAR-Fusion（https://zhuanlan.zhihu.com/p/586371542#h_586371542_6）

以第一个样本为例，5.fusion_from_fq/SRR11178348/

有一些文件和前面star比对结果类似就不再说，这里挑融合基因相关的看看

log文件

Log.final.out

                                  CHIMERIC READS:
                       Number of chimeric reads |       795505
                            % of chimeric reads |       2.06%

在STAR-Fusion的结果中，"CHIMERIC READS" 指的是被检测到具有融合特征的配对读取。融合基因事件通常由两个基因的某些部分融合形成，这可以导致在测序数据中检测到具有异常匹配的配对读取。这种异常匹配的配对读取就是 "CHIMERIC READS"，它表示两个不同的基因的片段在测序数据中被连在一起。这些 "CHIMERIC READS" 可能来自于基因融合点的跨越、基因重排等事件，是检测和识别融合基因的关键信息。根据这些异源的配对读取，STAR-Fusion可以进行融合事件的定位和推断。通常，"CHIMERIC READS" 数量的多少可以反映样本中融合基因事件的丰度或表达水平。具体分析和解读融合基因事件时，"CHIMERIC READS" 是一个重要的指标之一。

Chimeric.out.junction

我们这里看看和前面star调用相关参数得到的Chimeric.out.junction文件有什么区别，因为我们用前面的Chimeric.out.junction使用Kickstart模式的star-fusion并没有找到融合基因事件

$grep -v "^chr_" ../../4.mapping/SRR11178348Chimeric.out.junction | wc -l
379390
$grep -v "^chr_" Chimeric.out.junction | wc -l
1367384

可以发现我们现在Chimeric.out.junction文件包含了更多可能的融合基因事件记录

star-fusion.fusion_predictions.abridged.tsv

STAR-Fusion的输出是一个以制表符分隔的文件(star-fusion.fusion_predictions.tsv)和一个不包括证据融合reads识别的删减版本(star-fusion.fusion_predictions.abridged.tsv)

$head -2 star-fusion.fusion_predictions.abridged.tsv 
#FusionName JunctionReadCount SpanningFragCount est_J est_S SpliceType LeftGenLeftBreakpoint RightGene RightBreakpoint LargeAnchorSupport FFPM LeftBreakDinuc LeftBreakEntropy RightBreakDinuc RightBreakEntropy annots
E2F3--CDKAL1 48 9 48.00 9.00 ONLY_REF_SPLICE E2F3^ENSG00000112242.15 chr6:20488248:+CDKAL1^ENSG00000145996.11 chr6:20649293:+ YES_LDAS 1.4756 GT 1.8323 AG 1.9219 ["Klijn_CellLines","TCGA_StarF2019","CCLE_StarF2019","INTRACHROMOSOMAL[chr6:0.04Mb]","NEIGHBORS[40743]"]

2、获取上下游序列

下载基因组CDS、注释文件并格式化为外显子和CDS文件

Homo_sapiens.GRCh38.cds.all.fa.gz
Homo_sapiens.GRCh38.110.chr.gtf.gz

zcat Homo_sapiens.GRCh38.110.chr.gtf.gz |grep -w havana |perl -F"\t" -alne  '{next if $F[2] ne "exon";@a=split(/\"/,$F[8]); print join("\t",$a[1],$a[5],$a[9],$F[0],$F[3],$F[4],$F[4]-$F[3]+1) }' > Homo_sapiens.GRCh38.g2t2exon.txt

zcat Homo_sapiens.GRCh38.110.chr.gtf.gz |grep -w havana |perl -F"\t" -alne  '{next if $F[2] ne "CDS";@a=split(/\"/,$F[8]); print join("\t",$a[1],$a[5],$a[9],$F[0],$F[3],$F[4],$F[4]-$F[3]+1) }' > Homo_sapiens.GRCh38.g2t2cds.txt

关于gtf文件格式，参考：

文件格式 - GFF/GTF （https://www.bilibili.com/read/cv12909188）

grep -w havana：使用 grep 命令根据关键词 “havana” 进行匹配。-w 参数表示精确匹配整个单词。

$zcat Homo_sapiens.GRCh38.110.chr.gtf.gz | grep -v "^#" | awk '{print $2}' | sort | uniq -c
 122430 ensembl
 718375 ensembl_havana
    275 ensembl_havana_tagene
     21 ensembl_tagene
2404017 havana
 170147 havana_tagene
    144 insdc
   5637 mirbase

“havana” 的注释通常被认为是较为可靠和高质量的注释。因此，通过匹配 “havana” 可以过滤掉其他注释来源（如"ensembl"、"refseq"等），并仅保留具有 “havana” 注释来源的行。

打开文件仔细观察，可以发现CDS相比于exon序列，只是在transcript的5' 3'两端少了一部分

这里就涉及到一个非翻译区外显子的概念：

在鉴定lncRNA流程全套代码整理中，我们比较了CDS和mRNA的不同，主要区别在于CDS不包含两端的UTR部分

"exons and two codons"中的两个密码子当然就是起始密码子和终止密码子

那么CDS和exons片段序列严格意义上又有什么区别呢？在前面的文件里。可以发现CDS相比于exon序列，只是在transcript的5' 3'两端少了一部分，少的部分就是非翻译区exon

参考：原来外显子不是全都翻译啊！（https://zhuanlan.zhihu.com/p/87788625）

查看原始gtf文件也可进行区分，可以看看起始密码子位于哪个exon的哪个位置

后面提取序列进行可视化都以SRR11178348样本为例

将基因组坐标转为转录本坐标

library(stringr)

fusion <- read.table("/home/data/t120455/fusion_and_mut/5.fusion_from_fq/SRR11178348/star-fusion.fusion_pdictions.abridged.tsv")
head(fusion)
cds <- read.table("~/refdata/Homo_sapiens.GRCh38.g2t2exon.txt")
colnames(cds) <- c('gene','transcript','exon','chr','start','end','exon_length')
head(cds)

if (!exists("re")) {
  re <- array(dim = c(0, 9))  # 占位

if (!exists("re")) {
  re <- array(dim = c(0, 9))  # 占位
  colnames(re) <- c('gene','transcript','exon','chr','start','end','exon_length','pos','new_pos')
}
re

for (m in 1:nrow(fusion)) {
  x <- fusion[m,]
  g <- str_split(x[7],"\\^", simplify = TRUE)[1,2]
  g <- str_split(g,"\\.", simplify = TRUE)[1,1]
  pos <- as.numeric(str_split(x[8], ":", simplify = TRUE)[1,2])
  info <- cds[cds$gene == g,] # 有的基因有多个转录本
  if (nrow(info) > 0) {
    for (j in 1:length(split(info,info[,2]))) {
      y <- split(info,info[,2])[[j]]
      for (i in 1:nrow(y)) {
        z <- as.character(y[i,])
        if(z[5] <= pos & z[6] >= pos ){
          new <- sum(as.numeric(y[y[,3] < as.numeric(z[3]),7])) + pos - as.numeric(z[5])+1 
          re <- rbind(re,c(z,pos,new))
        }
        i <- 1 + 1
      }
      j <- j + 1
    }
  }
    m <- m + 1
}

head(re)

     gene              transcript        exon chr  start      end        exon_length pos        new_pos
[1,] "ENSG00000060237" "ENST00000535572" "25" "12" "907847"   "908034"   "188"       "908034"   "7060" 
[2,] "ENSG00000060237" "ENST00000537603" "2"  "12" "907847"   "908201"   "355"       "908034"   "386"  
[3,] "ENSG00000060237" "ENST00000540885" "1"  "12" "907720"   "908034"   "315"       "908034"   "315"  
[4,] "ENSG00000060237" "ENST00000675631" "24" "12" "907847"   "908034"   "188"       "908034"   "5879" 
[5,] "ENSG00000060237" "ENST00000676347" "18" "12" "907847"   "908034"   "188"       "908034"   "4278" 
[6,] "ENSG00000204138" "ENST00000463428" "2"  "1"  "28407410" "28407463" "54"        "28407463" "141"

根据转录本坐标获取转录本序列

fa <- readLines("~/refdata/Homo_sapiens.GRCh38.cds.all.fa.gz")
fa <- paste0(fa, collapse = "")
fa <- str_split(fa, ">")[[1]]
fa <- fa[-1]
tid <- str_split(fa, "\\.", simplify = TRUE)[,1]
head(tid)

seq <- str_split(fa, "]", simplify = TRUE)[,2]
head(seq)

re <- as.data.frame(re)
re <- re[re$transcript %in% tid,]
re$tseq <- seq[match(re$transcript,tid)]
head(re)
re$up <- unlist(lapply(1:nrow(re),function(i){
  #i=1
  new_pos <- as.numeric(re[i,9])
  l <- nchar(re[i,10])
  if(new_pos>500){
    return(substring(re[i,10],new_pos-500,new_pos))
  }else{
    return(substring(re[i,10],0,new_pos))
  }
  
}))
re$down <- unlist(lapply(1:nrow(re),function(i){
  #i=1
  new_pos <- as.numeric(re[i,9])
  l <- nchar(re[i,10])
  if((l-new_pos) >500){
    return(substring(re[i,10],new_pos,new_pos+500))
  }else{
    return(substring(re[i,10],new_pos,l))
  }
  
}))

要注意，有的转录本不存在于FASTA文件中，所以最后会损失一部分转录本无法获取序列

3、chimeraviz可视化

chimeraviz支持来自9种不同融合发现工具（deFuse、EricScript、InFusion、JAFFA、FusionCatcher、FusionMap、PRADA、SOAPfuse和STAR-FUSION）的输入

原推文是使用内置数据集来演示的，我们这里尝试能不能对我们上面自己找到的融合基因事件可视化

全局可视化

示例代码：


library(chimeraviz)

star_fusion_candidates <- system.file("extdata",
                                      "star-fusion.fusion_candidates.final.abridged.txt",
                                      package="chimeraviz")
fusions <- import_starfusion(star_fusion_candidates, "hg38")

length(fusions)
plot_circle(fusions)

> star_fusion_candidates
[1] "/home/data/t120455/R/x86_64-pc-linux-gnu-library/4.2/chimeraviz/extdata/star-fusion.fusion_candidates.final.abridged.txt"

打开发现和我们上面的得到的star-fusion.fusion_predictions.abridged.tsv格式一致

> getwd()
[1] "/home/data/t120455/fusion_and_mut/6.transcript_viz"

my_fusion_candidates <- c("../5.fusion_from_fq/SRR11178348/star-fusion.fusion_predictions.abridged.tsv")
my_fusions <- import_starfusion(my_fusion_candidates, "hg38")
length(my_fusions)
plot_circle(my_fusions)

单独可视化

示例代码：

需要融合位点附近的bam序列

原推文代码中存在部分代码不完整，且一些函数、参数需要更新，下面是我在原代码基础上修改后的代码和结果

defuse833ke <- system.file("extdata",
                           "defuse_833ke_results.filtered.tsv",
                           package = "chimeraviz")
fusions <- import_defuse(defuse833ke, "hg19", 1)
fusion <- get_fusion_by_id(fusions, 5267)
edbSqliteFile <- system.file("extdata",
                             "Homo_sapiens.GRCh37.74.sqlite",
                             package="chimeraviz")
edb <- ensembldb::EnsDb(edbSqliteFile)
                  
fusion5267and11759reads <- system.file("extdata",
                                       "fusion5267and11759reads.bam",
                                       package = "chimeraviz")
# 基于转录本
plot_fusion(fusion = fusion,
           bamfile = fusion5267and11759reads,
           edb = edb,
           non_ucsc = TRUE)

bamfile5267 <- system.file("extdata",
                           "5267readsAligned.bam",
                           package = "chimeraviz")

# 基于基因
plot_fusion(fusion = fusion,
           bamfile = bamfile5267,
           edb = edb,
           non_ucsc = TRUE,
           reduce_transcripts = TRUE)

可视化我们的结果第一个样本中第一个融合基因事件：

观察示例代码plot_fusion函数，有三个主要输入

fusion

这个最简单，我们前面已经通过import_starfusion读取了融合基因事件结果文件

my_fusion_candidates <- c("../5.fusion_from_fq/SRR11178348/star-fusion.fusion_predictions.abridged.tsv")
my_fusions <- import_starfusion(my_fusion_candidates, "hg38")

在示例代码中，使用了 get_fusion_by_id 通过id来获取单独的融合基因事件，我们这里就直接对列表操作[[1]]取第一个

> edbSqliteFile
[1] "/home/data/t120455/R/x86_64-pc-linux-gnu-library/4.2/chimeraviz/extdata/Homo_sapiens.GRCh37.74.sqlite"

edb <- ensembldb::EnsDb(edbSqliteFile)

> edb
EnsDb for Ensembl:
|Backend: SQLite
|Db type: EnsDb
|Type of Gene ID: Ensembl Gene ID
|Supporting package: ensembldb
|Db created by: ensembldb package from Bioconductor
|script_version: 0.0.1
|Creation time: Sat May 28 16:01:27 2016
|ensembl_version: 74
|ensembl_host: unknown
|Organism: Homo_sapiens
|genome_build: GRCh37
|DBSCHEMAVERSION: 1.0
|source_file: Homo_sapiens.GRCh37.74_subset.gtf
| No. of genes: 4.
| No. of transcripts: 21.

当有不认识的文件时，拿文件名去google一下是一个不错的选择

我这里检索时把37换成了38，因为我们自己使用的参考基因组是GRCh38

https://www.bioconductor.org/packages/devel/bioc/vignettes/ensembldb/inst/doc/ensembldb.html

生成和使用基于ENEMBL的注释软件包，看起来是使用ENSEMBL gtf文件的一个R包

简单看一下介绍：

ensembldb包提供了创建和使用以转录本为中心的注释数据库/包的功能。数据库的注释是使用 Perl API直接从 Ensembl获取的。除了从数据库中检索所有基因/转录模型和注释之外，该ensembldb 包还提供了一个过滤器框架，允许检索特定条目的注释，例如编码在基因上的基因lincRNA 基因的染色体区域或转录模型。该包的另一个主要目标是生成版本化注释包，即为特定 Ensembl 版本构建的注释包，并且也根据该版本命名（例如，用于 EnsDb.Hsapiens.v86Ensembl 代码数据库版本 86 的人类基因定义）。这确保了可重复性，因为如果有较新版本的注释包/版本可用，它也允许从特定 Ensembl 版本加载注释。它还允许同时加载多个注释包，以便比较 Ensembl 版本之间的基因模型。在下面的示例中，我们为智人加载基于 Ensembl 的注释包，Ensembl 版本 86。EnsDb提供对底层 SQLite 数据库的访问的对象绑定到变量 name EnsDb.Hsapiens.v86。

img

人GHCh38版本，就用这个了，我没对这个包的其他功能有什么探索，ENSEMBL的注释包我平时使用较多的是biomart

示例代码使用的数据集Homo_sapiens.GRCh37.74.sqlite只包含了4个genes，21个transcripts

EnsDb.Hsapiens.v86明显是汇总的

bamfile

需要融合位点附近的bam序列

关于sam/bam文件格式，参考：

SAM文件与Samtools （http://47.57.89.98/detail/85/）序列匹配中的CIGAR （http://47.57.89.98/detail/84/）

我这里直接拿star-fusio产生的Aligned.out.bam试试看

> plot_fusion(fusion = my_fusions[[1]],
+             bamfile = "../5.fusion_from_fq/SRR11178348/Aligned.out.bam",
+             edb = edb38,
+             non_ucsc = TRUE,
+             reduce_transcripts = TRUE)
Fetching transcripts for gene partners..
..transcripts fetched.
Fusion is intrachromosomal..
..but too far apart. Plot separate!
Selecting transcripts for E2F3..
..found transcripts of type exonBoundary
Selecting transcripts for CDKAL1..
..found transcripts of type exonBoundary
Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'x' in selecting a method for function 'colnames': error in evaluating the argument 'x' in selecting a method for function 'mcols': Unable to find index for BAM file '../5.fusion_from_fq/SRR11178348/Aligned.out.bam'. You can build an index using the following command:
 library(Rsamtools)
 indexBam("../5.fusion_from_fq/SRR11178348/Aligned.out.bam")

告诉我需要索引，建索引：

$samtools index Aligned.out.bam
[E::hts_idx_push] Chromosome blocks not continuous
[E::sam_index] Read 'SRR11178348.231796' with ref_name='chr2', ref_length=242193529, flags=99, pos=111858385 cannot be indexed
samtools index: failed to create index for "Aligned.out.bam"

报错解决，参考：

Samtools Index: Chromosome Blocks not Continuous （https://bioinformatics.stackexchange.com/questions/13997/samtools-index-chromosome-blocks-not-continuous）

于是先默认按照序列位置排序后再建索引：

$samtools sort Aligned.out.bam -o Aligned.out.sorted.bam
[bam_sort_core] merging from 5 files and 1 in-memory blocks...

$samtools index Aligned.out.sorted.bam

完整代码：

# install.packages("EnsDb.Hsapiens.v86_2.99.0.tar.gz",repos = NULL)
library(EnsDb.Hsapiens.v86)
## Making a "short cut"
edb38 <- EnsDb.Hsapiens.v86
## print some informations for this package
edb38
# 基于转录本
plot_fusion(fusion = my_fusions[[1]],
            bamfile = "../5.fusion_from_fq/SRR11178348/Aligned.out.sorted.bam",
            edb = edb38,
            non_ucsc = TRUE)
# 基于基因
plot_fusion(fusion = my_fusions[[1]],
            bamfile = "../5.fusion_from_fq/SRR11178348/Aligned.out.sorted.bam",
            edb = edb38,
            non_ucsc = TRUE,
            reduce_transcripts = TRUE)

参数说明：

reduce_transcripts
Boolean indicating whether or not to reduce all transcripts into a single transcript for each partner gene.

基于转录本

基于基因

我们这里使用的bam文件为star-fusion直接处理fq文件产生的

从大小上可以看出，tar-fusion直接处理fq文件产生的比对上的序列bam和star比对产生的bam文件不一样

sorted.bam相差了0.4G

这里也拿前面star产生的bam文件可视化看看：

基于转录本

基于基因

其实可以发现，对于我们目前这个E2F3-CDKAL1融合事件来说，并没有区别 57(48, 9)

参考： Bioconductor包chimeraviz嵌合RNA可视化（http://www.360doc.com/content/21/0714/12/76149697_986499751.shtml）顶部：显示融合的染色体位置。 支持断裂点（红色曲线）的discordant reads数10（其中split的6，spanning的4），注释的转录本及read数图。

简单看看这两个基因相关的报道

Medreading：

GeneCards：

E2F3得分最高的GeneHancer的Gene Targets中存在CDKAL1，可以为后面如果深入研究提供方向思路

gatk找突变

如果有读者仔细看过RNA-seq结题报告，就会发现在定量分析以外通常还会有SNP和INDEL分析。目前，对人类测序数据找突变最常用的软件是GATK，除了速度慢以外，没有其他明显缺点（可以通过部署Spark提高速度；当然，如果有钱，可以购买Sentieon，快了15-20倍）。和WES不同，RNA-seq对于外显子区域的覆盖度极度不均一，并且由于其数据来自转录本，跨越剪切位点的对比方式也使其与常规全外显子测序大不相同。对于RNA-seq与WES进行变异检测的异同，黄树嘉老师在《RNA-Seq能替代WES完成外显子的变异检测吗？》一文中有详细解释。本文后续内容主要来自GATK官网，但是这一部分的分析需要对WES检测变异具有一定的基础。如果跟不上下面的代码，墙裂推荐生信技能树的《肿瘤外显子专栏》，对于变异检测进行一个系统的了解之后再来follow。

这个专栏我曾经跟过：https://www.yuque.com/biotrainee/wes/mt7uba，可惜在不同 vcf 注释结果转 maf的时候发现注释结果有点问题，就没做下去

软件安装以及数据下载：

在语雀专栏中：

GATK，这里下载的版本为最新版本 4.1.4.1（截止时间 2019年 12月 31日），方法是：先进入 gatk 的官网 https://software.broadinstitute.org/gatk/，点击 Download (注：由于 GATK 官网界面的更新，新版本的界面已经不再如本文描述一样，且 ftp 链接也于 2020年6月1日关闭，所有数据库文件都迁移至谷歌云存储)

所以有些下载地址和方法存在变化，我这里主要参考我自己当时跟学时的笔记

作为一个菜鸟gatk我用的并不多，但每次使用都会阵痛，因为感觉对于一个小白来说入手不是很友好，有很多的依赖数据库文件、很多的参数，之前看到王通老师一篇推文 是时候将GATK彻底埋葬了 ，当然大佬看重的更多是软件性能和效果，我这等菜鸟还没到达那个层次

所以这里希望尽量整理清楚一些

官网下载安装

https://github.com/broadinstitute/gatk/releases

github release的文件一般是下载龟速，复制链接直接去服务器下更是，稍微好一点可以本地魔法下下来传给服务器，有时也比较慢

我这里检索加速下载的工具

参考：

下载速度慢？教你如何高速下载 github release！（https://zhuanlan.zhihu.com/p/521472595）

把文件地址输入进去，拿到加速链接用IDM下

传给服务器，解压打开，可以看到gatk可执行程序

$./gatk --version
Using GATK jar /home/data/t120455/gatk-4.3.0.0/gatk-package-4.3.0.0-local.jar
Running:
    java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /home/data/t120455/gatk-4.3.0.0/gatk-package-4.3.0.0-local.jar --version
The Genome Analysis Toolkit (GATK) v4.3.0.0
HTSJDK Version: 3.0.1
Picard Version: 2.27.5

conda安装

$mamba search gatk4 | tail -3
gatk4                        4.2.6.1  py36hdfd78af_1  anaconda/cloud/bioconda
gatk4                        4.3.0.0  py36hdfd78af_0  anaconda/cloud/bioconda
gatk4                        4.4.0.0  py36hdfd78af_0  anaconda/cloud/bioconda

$mamba search gatk | tail -3
gatk                             3.8          py36_2  anaconda/cloud/bioconda
gatk                             3.8          py36_3  anaconda/cloud/bioconda
gatk                             3.8          py36_4  anaconda/cloud/bioconda

这里可以看到gatk4.4对我们一些常用的包，如samtools和fastp版本进行了要求

Package       Version  Build           Channel                                                       Size
─────────────────────────────────────────────────────────────────────────────────────────────────────────────
  Install:
─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  + gatk4       4.4.0.0  py36hdfd78af_0  mirrors.bfsu.edu.cn/anaconda/cloud/bioconda                  364MB

  Upgrade:
─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  - lcms2          2.12  hddcbb42_0      mirrors.bfsu.edu.cn/anaconda/cloud/conda-forge                  
  + lcms2          2.14  h6ed2654_0      mirrors.bfsu.edu.cn/anaconda/cloud/conda-forge               262kB
  - libtiff       4.2.0  hf544144_3      mirrors.bfsu.edu.cn/anaconda/cloud/conda-forge                  
  + libtiff       4.4.0  h82bc61c_5      mirrors.bfsu.edu.cn/anaconda/cloud/conda-forge               485kB
  - openjdk      11.0.1  h516909a_1016   mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge         
  + openjdk      17.0.3  h58dac75_5      mirrors.bfsu.edu.cn/anaconda/cloud/conda-forge               169MB

  Downgrade:
─────────────────────────────────────────────────────────────────────────────────────────────────────────────

  - fastp        0.23.4  hadf994f_2      bioconda                                                        
  + fastp        0.23.3  h5f740d0_0      mirrors.bfsu.edu.cn/anaconda/cloud/bioconda                    3MB
  - htslib         1.17  h81da01d_2      mirrors.bfsu.edu.cn/anaconda/cloud/bioconda                     
  + htslib        1.3.1  h0cf4675_7      mirrors.bfsu.edu.cn/anaconda/cloud/bioconda                    1MB
  - libdeflate     1.18  h0b41bf4_0      mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge         
  + libdeflate     1.14  h166bdaf_0      mirrors.bfsu.edu.cn/anaconda/cloud/conda-forge                83kB
  - samtools       1.17  hd87286a_1      mirrors.bfsu.edu.cn/anaconda/cloud/bioconda                     
  + samtools        1.6  hc3601fc_10     mirrors.bfsu.edu.cn/anaconda/cloud/bioconda                 Cached

  Summary:

  Install: 1 packages
  Upgrade: 3 packages
  Downgrade: 4 packages

  Total download: 538MB

gatk4.3和我们目前的环境（stringtie鉴定lncRNA、star-fusion鉴定融合基因）是很兼容的

 + gatk4  4.3.0.0  py36hdfd78af_0  mirrors.bfsu.edu.cn/anaconda/cloud/bioconda     294MB

  Summary:

  Install: 1 packages

  Total download: 294MB

我们这里也选择下载gatk4.3，和前面官网下载的版本保持一致，随着一些软件和包的坑踩多了我也不喜欢装最新版

(lncRNA) t120455 00:42:14 ~ 
$which gatk
/home/data/t120455/miniconda3/envs/lncRNA/bin/gatk
(lncRNA) t120455 00:42:17 ~ 
$gatk --version
Using GATK jar /home/data/t120455/miniconda3/envs/lncRNA/share/gatk4-4.3.0.0-0/gatk-package-4.3.0.0-local.jar
Running:
    java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /home/data/t120455/miniconda3/envs/lncRNA/share/gatk4-4.3.0.0-0/gatk-package-4.3.0.0-local.jar --version
The Genome Analysis Toolkit (GATK) v4.3.0.0
HTSJDK Version: 3.0.1
Picard Version: 2.27.5

和官网安装的是一样的，其实我最希望conda安装帮我解决的是相关依赖数据库文件

数据库文件下载

我们这里先不着急往后走本项目gatk寻找突变

我们着重探索gatk数据库文件下载的最佳路径

所以这里我们根据语雀肿瘤外显子专栏测试下载提到的所需数据库文件：

https://www.yuque.com/biotrainee/wes/mt7uba

## gatk
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp_146.hg38.vcf.gz & 
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp_146.hg38.vcf.gz.tbi & 
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz & 
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.tbi & 
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.gz & 
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.fai & 
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.dict & 
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz & 
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi & 
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/funcotator/funcotator_dataSources.v1.6.20190124s.tar.gz &

这里给的是ftp server文件路径，不知道还能不能用

Resource bundle(https://gatk.broadinstitute.org/hc/en-us/articles/360035890811-Resource-bundle)

在官网这篇关于数据库文件下载的文章中提到：

FTP 服务器访问将很快被禁用，并且使用 FTP 文件路径的代码必须在 2020 年 6 月 1 日之前使用 Google Bucket 文件路径进行更新。

FTP 服务器适用于希望下载文件并在本地运行的人们。然而，FTP 是 Broad Institute 本地的（没有镜像），它对并发下载有严格的限制，并且一些国家的用户报告由于防火墙等原因难以访问它。由于这些（和其他）原因，FTP 服务器访问将于 2020 年 6 月 1 日之前被禁用。

此外还有两种assess the resource bundle的方法：Google Cloud bucket和Cromwell on Azure

下面以1000G_phase1.snps.high_confidence.hg38.vcf.gz文件（1.8G）为例进行下载测试，请注意三个方法并不是所有文件都包含，详细内容参考上面的Resource bundle(https://gatk.broadinstitute.org/hc/en-us/articles/360035890811-Resource-bundle)，

ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz

看起来还能用

wget:

$wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz
--2023-09-01 01:20:29--  ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz
           => ‘1000G_phase1.snps.high_confidence.hg38.vcf.gz’
Resolving ftp.broadinstitute.org (ftp.broadinstitute.org)... 69.173.70.223
Connecting to ftp.broadinstitute.org (ftp.broadinstitute.org)|69.173.70.223|:21... connected.
Logging in as gsapubftp-anonymous ... Logged in!
==> SYST ... done.    ==> PWD ... done.
==> TYPE I ... done.  ==> CWD (1) /bundle/hg38 ... done.
==> SIZE 1000G_phase1.snps.high_confidence.hg38.vcf.gz ... 1888262073
==> PASV ... done.    ==> RETR 1000G_phase1.snps.high_confidence.hg38.vcf.gz ... done.
Length: 1888262073 (1.8G) (unauthoritative)

hg38.vcf.gz                 6%[=>                                    ] 121.25M   335KB/s    eta 73m 11s

axel -n20:

速度还是比较客观

Google Buckets

https://console.cloud.google.com/storage/browser/genomics-public-data/resources/broad/hg38/v0

https://console.cloud.google.com/storage/browser/_details/genomics-public-data/resources/broad/hg38/v0/1000G_phase1.snps.high_confidence.hg38.vcf.gz

$wget -c https://console.cloud.google.com/storage/browser/_details/genomics-public-data/resources/broad/hg38/v0/1000G_phase1.snps.high_confidence.hg38.vcf.gz
--2023-09-01 01:30:00--  https://console.cloud.google.com/storage/browser/_details/genomics-public-data/resources/broad/hg38/v0/1000G_phase1.snps.high_confidence.hg38.vcf.gz
Resolving console.cloud.google.com (console.cloud.google.com)... 142.251.42.238
Connecting to console.cloud.google.com (console.cloud.google.com)|142.251.42.238|:443... connected.

$axel -n20 https://console.cloud.google.com/storage/browser/_details/genomics-public-data/resources/broad/hg38/v0/1000G_phase1.snps.high_confidence.hg38.vcf.gz
Initializing download: https://console.cloud.google.com/storage/browser/_details/genomics-public-data/resources/broad/hg38/v0/1000G_phase1.snps.high_confidence.hg38.vcf.gz

都连不上

gsutil：

之前我下载GTEx的数据库文件时用谷歌云推荐的gsutil，速度还行

下载gsutil

参考： gsutil 安装及使用（https://www.cnblogs.com/biostat-yu/p/15770784.html）

curl -O https://dl.google.com/dl/cloudsdk/channels/rapid/downloads/google-cloud-sdk-367.0.0-linux-x86_64.tar.gz
tar zxvf google-cloud-sdk-367.0.0-linux-x86_64.tar.gz

bash install.sh

根据提示信息，安装程序修改了我们的.bashrc配置文件

我们也把gsutil可执行程序目录路径放到环境变量中

$gsutil cp -r gs://genomics-public-data/resources/broad/hg38/v0/1000G_phase1.snps.high_confidence.hg38.vcf.gz ./
Copying gs://genomics-public-data/resources/broad/hg38/v0/1000G_phase1.snps.high_confidence.hg38.vcf.gz...
==> NOTE: You are downloading one or more large file(s), which would          
run significantly faster if you enabled sliced object downloads. This
feature is enabled by default but requires that compiled crcmod be
installed (see "gsutil help crcmod").

/ [0 files][409.4 MiB/  1.8 GiB]   24.2 MiB/s

需要注意这里gsutil下载必须给定输出路径

gsutil下载google buckets速度是axel -n20下载ftp速度的两倍

Azure

参考： 如何下载GATK resource bundle

下载安装azcopy，需要注意一些推文中的代码如果直接复制粘贴使用会报错，原因是格式问题，需要注意，这在我学习不同项目时遇到过

Azure：

https://learn.microsoft.com/en-us/azure/open-datasets/dataset-gatk-resource-bundle

这个页面很清爽，方便检索，还有文件的介绍信息

本质上是ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/的镜像，我之前用过1000G phase3做LD panel

把ftp换成https协议http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/就可以访问了

安装参考推文的方法将地址和token值拼接起来下载

对于1000G的文件，我这里测试West US 2的地址不行，West Central US可以

但是我没找到如何具体下载某个文件，比如1000G_phase1.snps.high_confidence.hg38.vcf.gz

介绍页面说“This dataset contains approximately 815 TB of data and is updated daily.”

有点害怕给服务器搞崩了，这里连上后下一会儿就给它取消了

简单看一下下载下来的目录

t120455 19:04:50 ~/refdata/dataset/phase1 
$tree
.
├── analysis_results
│   ├── integrated_call_sets
│   └── shapeit2_phased_haplotypes
│       └── ALL.chr9.SHAPEIT2_integrated_phase1_v3.20101123.snps_indels_svs.genotypes.all.vcf.gz.tbi
└── technical
    └── other_exome_alignments.alignment_indices

没找到如何具体下载某个文件，所以这个方法目前对我来说并不适用，希望知道的同学老师可以在评论区告诉我【Q1】

根据gatk找突变我们提前下载好所需的数据库文件，经过上面的测试，我们使用gsutil进行下载

但有时也很慢，这就具体情况具体选择方法了

需要注意的是dbsnp文件，参考语雀专栏评论区：

但是可能额因为vcf文件较大（超过10G），gsutil返回了这样的提示信息：

$gsutil cp -r gs://genomics-public-data/resources/broad/hg38/v0/Homo_sapiens_assembly38.dbsnp138.vcf ./
Copying gs://genomics-public-data/resources/broad/hg38/v0/Homo_sapiens_assembly38.dbsnp138.vcf...
==> NOTE: You are downloading one or more large file(s), which would          
run significantly faster if you enabled sliced object downloads. This
feature is enabled by default but requires that compiled crcmod be
installed (see "gsutil help crcmod").

CommandException: 
Downloading this composite object requires integrity checking with CRC32c,
but your crcmod installation isn't using the module's C extension, so the
hash computation will likely throttle download performance. For help
installing the extension, please see "gsutil help crcmod".

To download regardless of crcmod performance or to skip slow integrity
checks, see the "check_hashes" option in your boto config file.

NOTE: It is strongly recommended that you not disable integrity checks. Doing so
could allow data corruption to go undetected during uploading/downloading.

改用axel下载ftp

Downloaded 3.17687 Gigabyte(s) in 38:32 minute(s). (1440.56 KB/s)

interval_list

(lncRNA) t120455 00:50:30 ~/refdata/gatk_resource_bundle 
$gsutil cp -r gs://genomics-public-data/resources/broad/hg38/v0/wgs_calling_regions.hg38.interval_list ./

hapmap

gsutil cp -r gs://gcp-public-data--broad-references/hg38/v0/hapmap_3.3.hg38.vcf.gz ./

gsutil cp -r gs://gcp-public-data--broad-references/hg38/v0/1000G_omni2.5.hg38.vcf.gz ./

下面的处理我们仍以第一个样本SRR11178348为例进行处理，并使用conda下载下来的gatk4

区别于前面的工作目录命令脚本一般在第一层fusion_and_mut下，这里后面的脚本都统一放在了第二层7.gatk_mut下

1、向bam文件添加头文件并去重

gatk AddOrReplaceReadGroups \
-I ../4.mapping/SRR11178348Aligned.sortedByCoord.out.bam \
-O ./SRR11178348_rg_added_sorted.bam \
-ID SRR11178348 \
-LB rna \
-PL illumina \
-PU hiseq \
-SM SRR11178348 && \
gatk MarkDuplicates \
-I ./SRR11178348_rg_added_sorted.bam \
-O ./SRR11178348_dedup.bam  \
-M ./SRR11178348_dedup.metrics \
--CREATE_INDEX true \
--VALIDATION_STRINGENCY SILENT

2、转换MAPQ

gatk SplitNCigarReads \
-R ~/refdata/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa \
-I ./SRR11178348_dedup.bam \
-O ./SRR11178348_dedup_split.bam

A USER ERROR has occurred: Fasta index file file:///home/data/t120455/refdata/Homo_sapiens.GRCh38.dna_sm
.primary_assembly.fa.fai for reference file:///home/data/t120455/refdata/Homo_sapiens.GRCh38.dna_sm.prim
ary_assembly.fa does not exist. Please see http://gatkforums.broadinstitute.org/discussion/1601/how-can-
i-prepare-a-fasta-file-to-use-as-reference for help creating it.

参考： SAM文件与Samtools（http://47.57.89.98/detail/85/）

手动创建索引：

A USER ERROR has occurred: Fasta dict file file:///home/data/t120455/refdata/Homo_sapiens.GRCh38.dna_sm.
primary_assembly.dict for reference file:///home/data/t120455/refdata/Homo_sapiens.GRCh38.dna_sm.primary
_assembly.fa does not exist. Please see http://gatkforums.broadinstitute.org/discussion/1601/how-can-i-p
repare-a-fasta-file-to-use-as-reference for help creating it.

这个dict文件我在gatk的数据库文件里见过

.fasta文件的.dict文件通常是由Picard工具包中的CreateSequenceDictionary命令创建的，用于描述.fasta文件的参考序列

若直接改用gatk提供的参考基因组相关数据库文件，会出现这个报错：

 A USER ERROR has occurred: Input files reference and reads have incompatible contigs: No overlapping con tigs found.

reference contigs = [chr1, chr2, chr3, chr4, chr5, chr6, chr7, chr8, chr9, chr10, chr11, chr12, chr13, chr14, chr15, chr16, chr17, chr18, chr19, chr20, chr21, chr22, chrX, chrY, chrM, chr1_KI270706v1_random, chr1_KI270707v1_random, chr1_KI270708v1_random, chr1_KI270709v1_random, chr1_KI270710v1_random, chr1_KI270711v1_random, chr1_KI270712v1_random, chr1_KI270713v1_random, chr1_KI270714v1_random, chr2_KI270715v
1_random, chr2_KI270716v1_random, chr3_GL000221v1_random, ...]

 reads contigs = [1, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 2, 20, 21, 22, 3, 4, 5, 6, 7, 8, 9, MT, X, Y, KI270728.1, KI270727.1, KI270442.1, KI270729.1, GL000225.1, KI270743.1, ...]GL000008.2, GL000009.2, KI270747.1, KI270722.1, GL000194.1, KI270742.1, GL000205.2,

这个错误通常是因为参考序列文件和特征文件之间的染色体或序列名称不匹配导致的。当使用GATK等工具时，要求参考序列文件和特征文件（如vcf、bam文件）中的染色体/序列名称必须对应一致。

手动验证查看bam文件contig名称

$samtools view SRR11178348_dedup.bam | awk '{print $3}' | sort | uniq -c
6972121 1
2423838 10
3820170 11
4325551 12
1568207 13
3301455 14
2278217 15
3891309 16
4602996 17
 923565 18
4794034 19
5185623 2
2201475 20
3142953 21
1444028 22
3804606 3
2237300 4
2936031 5
4920276 6
4195418 7
2216511 8
2727988 9
   1944 GL000008.2
   3241 GL000009.2
   7321 GL000194.1
  14009 GL000195.1
   6639 GL000205.2
      2 GL000208.1
      6 GL000213.1
   5434 GL000214.1
     10 GL000216.2
   3730 GL000218.1
   6888 GL000219.1
1049361 GL000220.1
   2682 GL000221.1
   4577 GL000224.1
     88 GL000225.1
      2 KI270330.1
      2 KI270382.1
      2 KI270429.1
     88 KI270438.1
    134 KI270442.1
   7260 KI270466.1
  27229 KI270467.1
      2 KI270538.1
      8 KI270582.1
      2 KI270588.1
     20 KI270589.1
     34 KI270591.1
   2645 KI270706.1
     60 KI270707.1
    144 KI270708.1
     34 KI270709.1
  13411 KI270711.1
    118 KI270712.1
  10502 KI270713.1
     28 KI270714.1
     16 KI270715.1
    598 KI270716.1
   3784 KI270717.1
    176 KI270718.1
    152 KI270719.1
    420 KI270720.1
   1146 KI270721.1
    122 KI270722.1
     34 KI270725.1
     54 KI270726.1
     58 KI270727.1
    658 KI270728.1
      2 KI270729.1
    536 KI270731.1
      4 KI270732.1
1056468 KI270733.1
  10292 KI270734.1
      2 KI270735.1
      8 KI270736.1
   5126 KI270741.1
   4922 KI270742.1
    836 KI270743.1
   5388 KI270744.1
   2738 KI270745.1
    396 KI270746.1
    134 KI270747.1
     90 KI270748.1
     10 KI270749.1
     20 KI270750.1
   5306 KI270751.1
     64 KI270752.1
     16 KI270754.1
 373270 MT
1516002 X
 174341 Y

参考： Errors about input files having missing or incompatible contigs （https://gatk.broadinstitute.org/hc/en-us/articles/360035891131-Errors-about-input-files-having-missing-or-incompatible-contigs） GATK Input files reference and features have incompatible contigs: No overlapping contigs found. （https://www.biostars.org/p/9475152/）

在star比对的时候我们使用的是来自ENSEMBL的Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa

所以这里我们手动创建这个.dict文件，而gatk4相较于gatk已经整合了Picard工具包，所以这里我们直接用gatk4创建：

gatk CreateSequenceDictionary -R Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa -O Homo_sapiens.GRCh38.dna_sm.primary_assembly.dict

再次执行，结果文件：

gatk SplitNCigarReads \
-R ~/refdata/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa \
-I ./SRR11178348_dedup.bam \
-O ./SRR11178348_dedup_split.bam

img

3、进行碱基质量重校正

第一次尝试运行代码

gatk BaseRecalibrator \
-R ~/refdata/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa \
-I ./SRR11178348_dedup_split.bam \
--known-sites ~/refdata/gatk_resource_bundle/1000G_phase1.snps.high_confidence.hg38.vcf \
--known-sites ~/refdata/gatk_resource_bundle/dbsnp_146.hg38.vcf \
--known-sites ~/refdata/gatk_resource_bundle/Mills_and_1000G_gold_standard.indels.hg38.vcf \
-O ./SRR11178348_recal_data.table && \
gatk ApplyBQSR \
-R ~/refdata/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa \
-I ./SRR11178348_dedup_split.bam \
-bqsr ./SRR11178348_recal_data.table \
-O ./SRR11178348_BQSR.bam

又出现这个报错：

A USER ERROR has occurred: Input files reference and features have incompatible contigs: No overlapping contigs found.

reference contigs = [1, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 2, 20, 21, 22, 3, 4, 5, 6, 7, 8, 9, MT, X, Y, KI270728.1, KI270727.1, KI270442.1, KI270729.1, GL000225.1, KI270743.1, GL000008.2, GL000009.2, KI270747.1, KI270722.1, GL000194.1, KI270742.1, GL000205.2, GL000195.1, KI270736.1, KI270733.1, GL000224.1, GL000219.1, KI270719.1, GL000216.2, KI2707...]

features contigs = [chr1, chr2, chr3, chr4, chr5, chr6, chr7, chr8, chr9, chr10, chr11, chr12, chr13, chr14, chr15, chr16, chr17, chr18, chr19, chr20, chr21, chr22, chrX, chrY, chrM, chr1_KI270706v1_random, chr1_KI270707v1_random, chr1_KI270708v1_random, chr1_KI270709v1_random, chr1_KI270710v1_random, chr1_KI270711v1_random, chr1_KI270712v1_random, chr1_KI270713v1_random, chr1_KI

奇了怪了bam文件里的contigs名称一路走来怎么变了？并且和从gatk下载的fa一样，我最后用的参考基因组明明是一开始从ENSEMBL下载的啊

手动验证查看bam文件contig名称

$samtools view SRR11178348_dedup_split.bam | awk '{print $3}' | sort | uniq -c
9234929 1
3211308 10
5379127 11
6319267 12
2082983 13
4311425 14
3127047 15
5295168 16
6460004 17
1236139 18
7438580 19
6976103 2
2960933 20
3395504 21
1944173 22
5108098 3
2912255 4
3882250 5
6238225 6
5590652 7
2958436 8
3627543 9
   1951 GL000008.2
   3278 GL000009.2
   7339 GL000194.1
  14424 GL000195.1
   6775 GL000205.2
      2 GL000208.1
      6 GL000213.1
   5499 GL000214.1
     10 GL000216.2
   3761 GL000218.1
   7200 GL000219.1
1059917 GL000220.1
   2721 GL000221.1
   4602 GL000224.1
    107 GL000225.1
      2 KI270330.1
      2 KI270382.1
      2 KI270429.1
     88 KI270438.1
    145 KI270442.1
   7378 KI270466.1
  30200 KI270467.1
      2 KI270538.1
      8 KI270582.1
      2 KI270588.1
     20 KI270589.1
     42 KI270591.1
   2681 KI270706.1
     60 KI270707.1
    144 KI270708.1
     34 KI270709.1
  15846 KI270711.1
    118 KI270712.1
  10539 KI270713.1
     28 KI270714.1
     16 KI270715.1
    611 KI270716.1
   3784 KI270717.1
    176 KI270718.1
    152 KI270719.1
    472 KI270720.1
   1146 KI270721.1
    122 KI270722.1
     34 KI270725.1
     54 KI270726.1
     58 KI270727.1
    658 KI270728.1
      2 KI270729.1
    536 KI270731.1
      4 KI270732.1
1062111 KI270733.1
  10872 KI270734.1
      2 KI270735.1
      8 KI270736.1
   5203 KI270741.1
   5001 KI270742.1
    836 KI270743.1
   5449 KI270744.1
   2791 KI270745.1
    396 KI270746.1
    146 KI270747.1
     90 KI270748.1
     10 KI270749.1
     25 KI270750.1
   5395 KI270751.1
    112 KI270752.1
     16 KI270754.1
 373567 MT
2102382 X
 195505 Y

怎么和报错说的features contigs不一样？

哦，区别features于reads，这里的feature应该是指的是 --known-sites 参数指定的vcf文件，手动按照上面的方式查看contig name，证实想法，这也难怪，vcf文件都直接是gatk提供的

而今迈步从头越？我可是从star比对开始就用的从ENSEMBL下载的参考基因组，其实理论上一开始就用的gatk提供的参考基因组，后续使用gatk做其他分析就不会出现这些情况

能不能手动修改从gatk下载的数据库vcf文件的contigs名称？按理来说写脚本可以

问问chatgpt：

是的，你可以使用一些软件来直接修改VCF文件中的contig名称。以下是一些常用的工具和软件：

BCFtools：BCFtools是一个流行的工具集，可用于操作和处理VCF和BCF文件。使用 bcftools annotate命令可以修改VCF文件的contig名称。

bcftools annotate --rename-chrs <contig_mapping.txt> input.vcf -o output.vcf

其中，<contig_mapping.txt>是一个包含旧名称和新名称对应关系的文本文件。你需要提供一个符合格式要求的映射文件，将旧的contig名称映射到新的contig名称。

Picard工具包：Picard工具包中的 UpdateVcfSequenceDictionary命令可以用于更新VCF文件中的contig信息。

java -jar picard.jar UpdateVcfSequenceDictionary R=reference.fasta I=input.vcf O=output.vcf

在命令中，reference.fasta是与VCF文件关联的参考序列文件，input.vcf是输入的VCF文件，output.vcf是生成的包含更新contig信息的VCF文件。

VCFtools：VCFtools是一个广泛使用的VCF文件处理工具。你可以使用 vcf-annotate命令来修改VCF文件中的contig名称。

vcf-annotate --rename-chr <contig_mapping.txt> input.vcf > output.vcf

同样，你需要提供一个包含旧名称和新名称对应关系的映射文件 <contig_mapping.txt>。

请注意，在使用任何软件修改VCF文件之前，务必备份原始文件，以防出现意外情况。另外，确保映射文件中的名称对应正确，以避免错误的重命名。

发现有现成的工具不用自己写脚本

同时参考：

修改vcf文件中的染色体名（https://blog.csdn.net/liangbilin/article/details/108991510）

我们这里使用bcftools来进行修改：

mamba install bcftools

配置文件，根据实际情况修改，比如我这里chrMT改成chrM

修改名称

bcftools annotate --rename-chrs change_contig_name_vcf.config ./gatk_resource_bundle/1000G_phase1.snps.high_confidence.hg38.vcf -o ./gatk_resource_bundle/1000G_phase1.snps.high_confidence.hg38.RmChr.vcf

bcftools annotate --rename-chrs change_contig_name_vcf.config ./gatk_resource_bundle/Mills_and_1000G_gold_standard.indels.hg38.vcf -o ./gatk_resource_bundle/Mills_and_1000G_gold_standard.indels.hg38.RmChr.vcf

nohup bcftools annotate --rename-chrs change_contig_name_vcf.config ./gatk_resource_bundle/dbsnp_146.hg38.vcf -o ./gatk_resource_bundle/dbsnp_146.hg38.RmChr.vcf &

需要注意dbsnp_146.hg38.vcf

[W::vcf_parse] Contig 'chr1' is not defined in the header. (Quick workaround: index the file with tabix.)
Encountered an error, cannot proceed. Please check the error output above.
If feeling adventurous, use the --force option. (At your own risk!)

查看vcf文件可以看到确实没有像前面两个vcf文件那样的contig注释行

添加 --force 参数，也是同样的报错

参考：

生物信息问题杂烩（https://zhuanlan.zhihu.com/p/339566011?utm_id=0&wd=&eqid=daa8c63f0003925800000006648871ae#h_339566011_2）如果运行 bcftools 出现如下错误，说明需要对 vcf 创建索引文件 vcf文件处理&合并（https://zhuanlan.zhihu.com/p/650801709）最好不要用gzip直接压缩，可能会影响到后面建索引：我碰到的报错：index: "xxx.vcf.gz" is in a format that cannot be usefully indexed 所以压缩vcf文件最好用bgzip

创建索引后修改成功

bgzip dbsnp_146.hg38.vcf
bcftools index dbsnp_146.hg38.vcf.gz
cd ..
nohup bcftools annotate --rename-chrs change_contig_name_vcf.config ./gatk_resource_bundle/dbsnp_146.hg38.vcf.gz -o ./gatk_resource_bundle/dbsnp_146.hg38.RmChr.vcf.gz &
cd gatk_resource_bundle
gunzip dbsnp_146.hg38.RmChr.vcf.gz

我猜前面两个vcf文件因为有contig注释行，所以bcftools修改contig名称不需要索引

第二次尝试运行代码

gatk BaseRecalibrator \
-R ~/refdata/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa \
-I ./SRR11178348_dedup_split.bam \
--known-sites ~/refdata/gatk_resource_bundle/1000G_phase1.snps.high_confidence.hg38.RmChr.vcf \
--known-sites ~/refdata/gatk_resource_bundle/dbsnp_146.hg38.RmChr.vcf \
--known-sites ~/refdata/gatk_resource_bundle/Mills_and_1000G_gold_standard.indels.hg38.RmChr.vcf \
-O ./SRR11178348_recal_data.table && \
gatk ApplyBQSR \
-R ~/refdata/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa \
-I ./SRR11178348_dedup_split.bam \
-bqsr ./SRR11178348_recal_data.table \
-O ./SRR11178348_BQSR.bam

A USER ERROR has occurred: Input /home/data/t120455/refdata/gatk_resource_bundle/1000G_phase1.snps.high_
confidence.hg38.RmChr.vcf must support random access to enable queries by interval. If it's a file, please index it using the bundled tool IndexFeatureFile

A USER ERROR has occurred: Input /home/data/t120455/refdata/gatk_resource_bundle/dbsnp_146.hg38.RmChr.vcf must support random access to enable queries by interval. If it's a file, please index it using the bu
ndled tool IndexFeatureFile

A USER ERROR has occurred: Input /home/data/t120455/refdata/gatk_resource_bundle/Mills_and_1000G_gold_st
andard.indels.hg38.RmChr.vcf must support random access to enable queries by interval. If it's a file, p
lease index it using the bundled tool IndexFeatureFile

参考：

2021-03-31 为VCF文件建立索引（.idx) （https://www.jianshu.com/p/307bb7b83c3f） gatk IndexFeatureFile -I XXX.vcf

前面我们用bcftools给vcf.gz建立索引，而gatk可以直接对vcf建立索引

gatk IndexFeatureFile -I 1000G_phase1.snps.high_confidence.hg38.RmChr.vcf
gatk IndexFeatureFile -I dbsnp_146.hg38.RmChr.vcf

gatk IndexFeatureFile -I Mills_and_1000G_gold_standard.indels.hg38.RmChr.vcf

需要注意Mills_and_1000G_gold_standard.indels.hg38.RmChr.vcf

A USER ERROR has occurred: Error while trying to create index for Mills_and_1000G_gold_standard.indels.hg38.RmChr.vcf. Error was: htsjdk.tribble.TribbleException: The provided VCF file is malformed at approximately line number 3475: The reference allele cannot be missing

报错信息说第3475行的变异位点没有ref allele

查看：

t120455 23:46:56 ~/refdata/gatk_resource_bundle 
$sed -n "3475p" Mills_and_1000G_gold_standard.indels.hg38.RmChr.vcf 
1 0 . . . . . .
t120455 23:47:19 ~/refdata/gatk_resource_bundle 
$sed -n "3474,3476p" Mills_and_1000G_gold_standard.indels.hg38.RmChr.vcf 
1 965337 rs59448239 CTTAT C 37990.2 PASS set=Intersect1000GMinusSI
1 0 . . . . . .
t120455 23:47:33 ~/refdata/gatk_resource_bundle 
$sed -n "3474,3477p" Mills_and_1000G_gold_standard.indels.hg38.RmChr.vcf 
1 965337 rs59448239 CTTAT C 37990.2 PASS set=Intersect1000GMinusSI
1 0 . . . . . .
t120455 23:47:52 ~/refdata/gatk_resource_bundle 
$wc -l Mills_and_1000G_gold_standard.indels.hg38.RmChr.vcf
3475 Mills_and_1000G_gold_standard.indels.hg38.RmChr.vcf

可以发现就是最后一行，干脆删掉

t120455 23:48:08 ~/refdata/gatk_resource_bundle 
$sed -n "3475p" Mills_and_1000G_gold_standard.indels.hg38.RmChr.vcf 
1 0 . . . . . .
t120455 23:51:56 ~/refdata/gatk_resource_bundle 
$sed -i "3475d" Mills_and_1000G_gold_standard.indels.hg38.RmChr.vcf 
t120455 23:52:14 ~/refdata/gatk_resource_bundle 
$sed -n "3475p" Mills_and_1000G_gold_standard.indels.hg38.RmChr.vcf 
t120455 23:52:17 ~/refdata/gatk_resource_bundle 
$wc -l Mills_and_1000G_gold_standard.indels.hg38.RmChr.vcf
3474 Mills_and_1000G_gold_standard.indels.hg38.RmChr.vcf

创建索引成功

第三次运行代码

gatk BaseRecalibrator \
-R ~/refdata/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa \
-I ./SRR11178348_dedup_split.bam \
--known-sites ~/refdata/gatk_resource_bundle/1000G_phase1.snps.high_confidence.hg38.RmChr.vcf \
--known-sites ~/refdata/gatk_resource_bundle/dbsnp_146.hg38.RmChr.vcf \
--known-sites ~/refdata/gatk_resource_bundle/Mills_and_1000G_gold_standard.indels.hg38.RmChr.vcf \
-O ./SRR11178348_recal_data.table && \
gatk ApplyBQSR \
-R ~/refdata/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa \
-I ./SRR11178348_dedup_split.bam \
-bqsr ./SRR11178348_recal_data.table \
-O ./SRR11178348_BQSR.bam

成功，结果文件：

4、call mutations

gatk HaplotypeCaller \
-R ~/refdata/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa \
-I ./SRR11178348_BQSR.bam \
-O ./SRR11178348.vcf \
-L ~/refdata/gatk_resource_bundle/wgs_calling_regions.hg38.RmChr.interval_list

这里需要注意 -L 参数使用的interval_list文件

在GATK的HaplotypeCaller工具中，-L参数用于指定用于分析的区域列表文件路径。该参数允许指定特定的基因组区域，只对这些区域进行变异分析，而不是对整个基因组进行分析。这在特定感兴趣的区域或需要加快分析速度时非常有用，通过使用-L参数并提供区域列表文件路径，HaplotypeCaller工具将仅分析指定的区域，减少了分析的范围，从而节省了时间和计算资源。

参考：

Broad hg38 ICE interval list （https://gatk.broadinstitute.org/hc/en-us/community/posts/360057855032-Broad-hg38-ICE-interval-list） Intervals and interval lists （https://gatk.broadinstitute.org/hc/en-us/articles/360035531852-Intervals-and-interval-lists） How to choose interval_list? （https://gatk.broadinstitute.org/hc/en-us/community/posts/8622509405339-How-to-choose-interval-list-） You should use the wgs_calling_regions.hg38.interval_list for single-sample and joint calling.

文件下载：

我们前面使用的google bucket

https://console.cloud.google.com/storage/browser/genomics-public-data/resources/broad/hg38/v0

而在

How to choose interval_list? （https://gatk.broadinstitute.org/hc/en-us/community/posts/8622509405339-How-to-choose-interval-list- 和 Resource bundle(https://gatk.broadinstitute.org/hc/en-us/articles/360035890811-Resource-bundle)

中google bucket地址 https://console.cloud.google.com/storage/browser/gcp-public-data--broad-references/hg38/v0

看起来维护时间更新，数据库文件更多

下载interval_list

(lncRNA) t120455 00:50:30 ~/refdata/gatk_resource_bundle 
$gsutil cp -r gs://genomics-public-data/resources/broad/hg38/v0/wgs_calling_regions.hg38.interval_list ./

直接去掉chr

sed 's/chr//g'

运行代码，结果文件：

img

5、对mutations进行过滤

像前面一样，对vcf文件修改contig名称并构建索引

bcftools annotate --rename-chrs change_contig_name_vcf.config input.vcf -o input.RmChr.vcf

gatk IndexFeatureFile -I XX.RmChr.vcf

尝试运行代码：

model:SNP

gatk VariantRecalibrator \
-R ~/refdata/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa \
-V ./SRR11178348.vcf \
-resource:hapmap,known=false,training=true,truth=true,prior=15.0 ~/refdata/gatk_resource_bundle/hapmap_3.3.hg38.RmChr.vcf \
-resource:omini,known=false,training=true,truth=false,prior=12.0 ~/refdata/gatk_resource_bundle/1000G_omni2.5.hg38.RmChr.vcf \
-resource:1000G,known=false,training=true,truth=false,prior=10.0 ~/refdata/gatk_resource_bundle/1000G_phase1.snps.high_confidence.hg38.RmChr.vcf \
-resource:dbsnp,known=true,training=false,truth=false,prior=6.0 ~/refdata/gatk_resource_bundle/dbsnp_146.hg38.RmChr.vcf \
-an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum \
-mode SNP \
-tranche 100.0 -tranche 99.9 -tranche 99.0 -tranche 95.0 -tranche 90.0 \
--tranches-file SRR11178348.snps.tranches \
-O SRR11178348.snps.recal

gatk ApplyVQSR \
-R ~/refdata/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa \
-V ./SRR11178348.vcf \
--tranches-file SRR11178348.snps.tranches \
--recal-file SRR11178348.snps.recal \
-ts-filter-level 99.0 \
-mode SNP \
-O SRR11178348.snps.VQSR.vcf.gz

报错

13:40:16.601 INFO  ProgressMeter - Traversal complete. Processed 241459 total variants in 4.8 minutes.
13:40:16.641 INFO  VariantDataManager - DP:      mean = 7.30     standard deviation = 28.00
13:40:16.674 INFO  VariantDataManager - QD:      mean = 24.32    standard deviation = 7.13
13:40:16.691 INFO  VariantDataManager - FS:      mean = 0.25     standard deviation = 1.26
13:40:16.719 INFO  VariantDataManager - SOR:     mean = 1.31     standard deviation = 0.82
13:40:16.759 INFO  VariantDataManager - ReadPosRankSum:          mean = 0.02     standard deviation = 0.
96
13:40:16.787 INFO  VariantDataManager - MQRankSum:       mean = 0.00     standard deviation = 0.00
13:40:16.807 INFO  VariantRecalibrator - Shutting down engine
[September 3, 2023 1:40:16 PM CST] org.broadinstitute.hellbender.tools.walkers.vqsr.VariantRecalibrator 
done. Elapsed time: 4.82 minutes.
Runtime.totalMemory()=8761901056
***********************************************************************

A USER ERROR has occurred: Bad input: Found annotations with zero variance. They must be excluded before
 proceeding.

可以发现

13:40:16.787 INFO  VariantDataManager - MQRankSum:       mean = 0.00     standard deviation = 0.00

报错信息让我们排除掉这个注释

在GATK的VariantRecalibrator工具中，-an参数用于指定用于校准的注释特征。这些特征是基因组变异数据中的一些统计信息，用于帮助区分真实变异和假阳性。

以下是命令中使用的六个特征的解释：

DP（Read Depth）：该特征表示在某个位点上的总测序深度，即所有覆盖该位点的读取数量的总和。
QD（Quality by Depth）：该特征是变异质量（QUAL）除以测序深度（DP），用于衡量变异质量相对于测序深度的比值。较高的QD值通常表示更可靠的变异。
FS（Fisher Strand Bias）：该特征衡量单个位点附近的两种不同链的碱基的比例，以捕捉可能的测序偏差情况。较低的FS值表示较少的偏差。
SOR（Strand Odds Ratio）：该特征比较在正链和负链上观察到的变异读数的比例，以探测潜在的偏性。较低的SOR值表示较少的偏性。
ReadPosRankSum：该特征比较在变异位置上校正的碱基相对于变异位置外校正的碱基的质量得分。用于衡量碱基在读取中的位置对变异的贡献。较高的ReadPosRankSum值表示较高的质量。
MQRankSum（Mapping Quality Rank Sum Test）：该特征比较与该位点关联的变异碱基与参考碱基的比例，以衡量测序质量的偏差。较高的MQRankSum值表示较高的质量。

这些特征提供了关于基因组变异的不同方面的信息，通过将这些特征包含在质量校准中，可以更好地筛选和区分真实变异和假阳性。

参考：

GATK报错信息汇总（https://www.jianshu.com/p/dd60506b06a2）所使用的一些tags也就是annotations是不适用的，应该把这些注释去掉

去除后运行成功

mode:INDEL

gatk VariantRecalibrator \
-R ~/refdata/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa \
-V ./SRR11178348.snps.VQSR.vcf.gz \
-resource:mills,known=true,training=true,truth=true,prior=12.0 ~/refdata/gatk_resource_bundle/Mills_and_
1000G_gold_standard.indels.hg38.RmChr.vcf \
-an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum \
-mode INDEL \
--max-gaussians 6 \
--tranches-file SRR11178348.indels.tranches \
-O SRR11178348.snps.indels.recal

gatk ApplyVQSR \
-R ~/refdata/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa \
-V ./SRR11178348.snps.VQSR.vcf.gz \
--tranches-file SRR11178348.indels.tranches \
--recal-file SRR11178348.snps.indels.recal \
-ts-filter-level 99.0 \
-mode INDEL \
-O SRR11178348.VQSR.vcf.gz

报错

20:25:41.178 INFO  ProgressMeter - Traversal complete. Processed 241459 total variants in 0.0 minutes.
20:25:41.211 INFO  VariantDataManager - DP:   mean = 2.00  standard deviation = 0.00
20:25:41.242 INFO  VariantDataManager - QD:   mean = 27.73  standard deviation = 0.00
20:25:41.262 INFO  VariantDataManager - FS:   mean = -0.00  standard deviation = 0.00
20:25:41.274 INFO  VariantDataManager - SOR:   mean = 2.30  standard deviation = 0.00
20:25:41.299 INFO  VariantRecalibrator - Shutting down engine
[September 3, 2023 8:25:41 PM CST] org.broadinstitute.hellbender.tools.walkers.vqsr.VariantRecalibrator done. Elapsed time: 0.04 minutes.
Runtime.totalMemory()=2454716416
***********************************************************************

A USER ERROR has occurred: Bad input: Values for ReadPosRankSum annotation not detected for ANY training variant in the input callset. VariantAnnotator may be used to add these annotations.

A USER ERROR has occurred: Bad input: Values for MQRankSum annotation not detected for ANY training variant in the input callset. VariantAnnotator may be used to add these annotations.

去除这两个注释后，仍报错

20:45:32.352 INFO  ProgressMeter - Traversal complete. Processed 241459 total variants in 0.0 minutes.
20:45:32.370 INFO  VariantDataManager - DP:   mean = 2.00  standard deviation = 0.00
20:45:32.387 INFO  VariantDataManager - QD:   mean = 27.73  standard deviation = 0.00
20:45:32.395 INFO  VariantDataManager - FS:   mean = -0.00  standard deviation = 0.00
20:45:32.400 INFO  VariantDataManager - SOR:   mean = 2.30  standard deviation = 0.00
20:45:32.407 INFO  VariantRecalibrator - Shutting down engine
[September 3, 2023 8:45:32 PM CST] org.broadinstitute.hellbender.tools.walkers.vqsr.VariantRecalibrator done. Elapsed time: 0.05 minutes.
Runtime.totalMemory()=2667577344
***********************************************************************

A USER ERROR has occurred: Bad input: Found annotations with zero variance. They must be excluded before proceeding.

可以发现这四个注释的variance都为0

如果某个注释特征的方差为0，则意味着该特征在所有的样本中的取值都是相同的，没有变异性

为什么会出现这个情况？【Q2】

小结

以上就是本周全部内容，又是一个新的流程，后面可能没有太多时间做这样的事，因为debug太耗时间了，可以发现我上面基本上一直在debug，导致我目前还不能很好地解读代码参数和结果文件背后的生物学意义，这才是最重要的

这套流程可以和肿瘤外显子数据分析指南语雀专栏（https://www.yuque.com/biotrainee/wes/ipgs5k#z7xmO）比较学习

对于gatk的vcf结果文件，我们可以通过不同软件进行注释或转换成maf文件

原推文提到：

其实，在获得VCF文件之后，还可以使用数据库文件或自有的测序文件进行germline突变的去除，以进一步提高结果的可信度。

语雀专栏走的就是这个流程：

单样本找 Germline SNPS + INDELs 我们通常说的突变，可以分为种系突变 Germline mutations（又称生殖突变）和体细胞突变 Somatic mutations。前者是指在生殖细胞中发生的任何可检测、可遗传的突变。在生殖细胞系以外的细胞中发生的突变称为体细胞突变，也称作获得性突变，是不可遗传的。使用 GATK 的流程进行 Germline mutations，主要用到的工具是 HaplotypeCaller，流程示意图如下（最后得到的 vcf 文件应该还有走 VQSR 进行过滤） mutect 找 Somatic mutations 前面单样本走 GATK HaplotypeCaller 流程得到 vcf 文件主要包含了 germline mutations 信息。但是我们这里分析的是肿瘤样本，包含成对的 tumor 和 normal 样本，更加关心的是 tumor 样本中所包含的 somatic mutations 信息。所以我们应该走 GATK 的 somatic mutations 流程

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2023-09-05，如有侵权请联系 cloudcommunity@tencent.com 删除

数据库