一、外显子测序分析可以得到以下结果：

变异位点的检测和注释
基因突变的检测和注释
拼接剪接位点的检测和注释
基因外显子区域的覆盖度分析
基因外显子区域的差异表达分析
基因功能富集分析

二、利用外显子测序数据得到肿瘤特异性靶标的流程一般包括以下几个步骤：

数据预处理：包括测序数据的质量控制、去除低质量序列、比对到基因组参考序列等。
变异检测：利用不同的软件对测序数据进行变异检测，筛选出与肿瘤相关的变异位点。
变异注释：对变异位点进行注释，包括变异的类型、位置、影响等信息。
功能分析：对注释的结果进行功能富集分析，筛选出与肿瘤相关的基因和通路。
靶标筛选：根据功能分析的结果，筛选出与肿瘤相关的靶标，进行进一步的实验验证。

具体流程还需要根据具体的实验设计和数据分析要求进行调整。

外显子测序数据分析中常用的软件和工具：

数据预处理：Trimmomatic、FastQC、BWA、SAMtools等。
变异检测：GATK、VarScan、Mutect2、Strelka等。
变异注释：ANNOVAR、VEP、SNPEff等。
功能分析：DAVID、GOseq、KEGG等。
靶标筛选：Cytoscape、STRING、GeneMANIA等。

需要注意的是，具体软件的选择还要根据实验设计、数据类型和分析需求等因素进行综合考虑。

基础命令：
bwa mem hg38.fa HRR573094_f1.fq.gz HRR573094_r2.fq.gz > mem-pe.3094.sam

可添加的参数：

-t：指定线程数。可以根据计算机配置来设置。比如，有8个CPU核心可用，可以设置为"-t 8"，以加快比对速度。
-M：生成CIGAR字符串中的M标记，表示匹配和误配。这个选项会让BWA生成一个SAM文件，其中包含了所有的比对结果。如果只需要最好的比对结果，可以去掉这个选项。
-R：指定read group信息。这个信息对于后续的数据分析很重要，建议填写。可以参考以下格式：-R "@RG\tID:group1\tSM:sample1\tPL:illumina\tPU:unit1"
-B：指定比对算法。BWA有三种比对算法可以选择，分别是"mem"、"bwasw"和"aln"。其中，默认的是"mem"算法，适用于短读长于70bp的情况。

综上所述，可以使用以下命令进行比对：
bwa mem -t 8 -M -R "@RG\tID:group1\tSM:sample1\tPL:illumina\tPU:unit1" ref.fa read1.fq.gz read2.fq.gz > mem-pe.sam

使用循环语句来批量处理数据。下面是一个示例脚本，可以批量处理一个目录下的所有双端测序数据：

bash
#!/bin/bash
# 批量比对双端测序数据

# 1. 定义变量
REF="hg38.fa"
THREADS=8
DIR="data"  # 存放数据的目录
OUTDIR="output"  # 存放比对结果的目录

# 2. 创建输出目录
mkdir -p $OUTDIR

# 3. 处理每个样本
for file in $DIR/*.fq.gz
do

# 3.1 获取文件名和路径
filename=$(basename "$file")
sample="${filename%.*}"
   
 # 3.2 执行比对命令
bwa mem -t $THREADS -M -R "@RG\tID:$sample\tSM:$sample\tPL:illumina\tPU:unit1" $REF $DIR/${sample}_1.fq.gz $DIR/${sample}_2.fq.gz > $OUTDIR/${sample}.sam

# 3.3 转换为BAM格式并排序
samtools view -b -S $OUTDIR/${sample}.sam | samtools sort -o $OUTDIR/${sample}.bam -

# 3.4 建立索引
samtools index $OUTDIR/${sample}.bam
done

上述脚本中，我们首先定义了参考基因组文件和线程数等变量，然后创建了一个用于存放比对结果的目录。接下来使用for循环遍历指定目录下的所有双端测序数据，依次执行比对命令、转换为BAM格式并排序、建立索引等操作。

可以将上述脚本保存为.sh文件，然后使用终端运行即可。需要安装BWA和SAMtools等软件，并将它们的可执行文件添加到系统环境变量中，才能在终端中直接运行命令。

在Illumina测序中，每个测序文库都有一个唯一的标识符，称为“read group”。这个标识符可以用来区分不同的测序文库，并且可以提供有关测序数据的更多信息，例如测序平台、测序日期、测序文库类型等。在数据分析中，read group信息可以帮助我们更好地理解和解释数据。

下面是一个read group信息的示例：

@RG\tID:group1\tSM:sample1\tPL:illumina\tPU:unit1

其中，\t表示制表符，\n表示换行符。各个字段的含义如下：

ID：测序文库的唯一标识符。
SM：样本名称，即测序文库对应的样本名称。
PL：测序平台，例如Illumina、Ion Torrent等。
PU：测序文库的物理标识符，例如测序文库的barcode或lane信息等。

这个示例中的read group信息表明这个测序文库的唯一标识符是group1，对应的样本名称是sample1，测序平台是Illumina，物理标识符是unit1。

在Cell Ranger中，使用-R参数可以指定read group信息。如果您不确定如何填写read group信息，可以参考Illumina官方文档或者咨询测序服务商。

shell的命令搞不定。。自己写的python脚本挺好用。。

用bwa比对

#批量执行bwa
def bwa():
    import os
    for i in range(573200,573230):
        x = "HRR" + str(i)
        cmd_string = "bwa mem -t 8 -M hg38.fa "+x+"_f1.fq.gz "+x+"_r2.fq.gz > /mnt/bwa-0.7.17/output/"+x+".sam"
        print('x：{}'.format(cmd_string))
        print(os.popen(cmd_string).read())
                
bwa()

sam 转bam

conda install -c bioconda samtools

#批量执行samtools
def samtools():
    import os
    for i in range(573185,573200):
        x = "HRR" + str(i)
        cmd_string = "samtools view -bS /mnt/bwa-0.7.17/output/"+x+".sam > /mnt/bwa-0.7.17/output/"+x+".bam"
        print('x：{}'.format(cmd_string))
        print(os.popen(cmd_string).read())
samtools()

sort的命令终于能用了，网上好多教程是错的

排序

def sort():
    import os
    for i in range(573097, 573099):
        x = "HRR" + str(i)
        cmd_string = "samtools sort -@ 8 -m 8G -O bam -o /mnt/bwa-0.7.17/output/" + x + ".bam /mnt/bwa-0.7.17/output/" + x + ".bam"
        print('x：{}'.format(cmd_string))
        print(os.popen(cmd_string).read())
sort()

本来的报错：File "samtools"， line 1SyntaxError: Non-UTF-8 code starting with '\xde' in file samtools on line 2, but no encoding declared; see http://python.org/dev/peps/pep-0263/ for details我在第一行加了在coding: utf-8 --没有用，后来把能运行的脚本复制然后放这个的代码。

标记重复序列-去重

#批量执行sambamba
def bwa():
    import os
    for i in range(573098,573101):
        x = "HRR" + str(i)
        cmd_string = "sambamba markdup -t 8 /mnt/bwa-0.7.17/output/"+x+".bam /mnt/bwa-0.7.17/output/"+x+".marked.bam"
        print('x：{}'.format(cmd_string))
        print(os.popen(cmd_string).read())
bwa()

创建索引

进行sambamba的markdup，就会生成.bai文件，所以我想是不是已经有index了

index区域重比对，重比对区域校正

发现网上的教程也是五花八门各不相同。。所以就一个个试哪个能用

生成fai文件：
samtools faidx Homo_sapiens_assembly38.fasta

太久没wes了其实工作路径在/mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0

indel区域重比对

gatk BaseRecalibrator -I /mnt/bwa-0.7.17/output/HRR573099.marked.RG.bam -R /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/Homo_sapiens_assembly38.fasta --known-sites /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz --known-sites /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/dbsnp_146.hg38.vcf.gz --known-sites /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz -O /mnt/bwa-0.7.17/output/HRR573099.recal.table

出现报错“A USER ERROR has occurred: Number of read groups must be >= 1, but is 0”原因是没有read group分组信息。
教程：https://uteric.github.io/WES/wes01/
理应在bwa的时候就添加分组信息去跑（bwa mem -t 5 -R "@RG\tID:$sample\tSM:$sample\tLB:WGS\tPL:Illumina" $INDEX $fq1 $fq2 | samtools sort -@ 5 -o $sample.bam -），但是没有添加，可以用gatk AddOrReplaceReadGroups补救一下，示例：
gatk AddOrReplaceReadGroups -I /mnt/d/wes/bam/clean.bam -O /mnt/d/wes/bam/clean.RG.bam -LB genome -PL ILLUMINA -PU unit1 -SM sample1
尝试gatk AddOrReplaceReadGroups -I /mnt/bwa-0.7.17/output/HRR573099.marked.bam -O /mnt/bwa-0.7.17/output/HRR573099.marked.RG.bam -LB WGS -PL Illumina -PU unit1 -SM HRR573099成功，输入samtools view -H /mnt/bwa-0.7.17/output/HRR573099.marked.RG.bam | grep '@RG'反馈是：
@RG ID:1 LB:WGS PL:Illumina SM:HRR573099 PU:unit1

用循环运行：
/mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/gatk.py

def gatk():
    import os
    for i in range(573100,573131):
        x = "HRR" + str(i)
        cmd_string = "gatk AddOrReplaceReadGroups -I /mnt/bwa-0.7.17/output/"+x+".marked.bam -O /mnt/bwa-0.7.17/output/"+x+".marked.RG.bam -LB WGS -PL Illumina -PU unit1 -SM "+x
        print('x：{}'.format(cmd_string))
        print(os.popen(cmd_string).read())
                
gatk()

PCR去重

gatk MarkDuplicates -I /mnt/bwa-0.7.17/output/HRR573099.marked.RG.bam -O /mnt/bwa-0.7.17/output/HRR573099.nodup.bam -M /mnt/bwa-0.7.17/output/HRR573099.metrics.txt --REMOVE_DUPLICATES true --ASSUME_SORTED true
会生成HRR573099.nodup.bam、HRR573099.metrics.txt
批量：

def gatk():
    import os
    for i in range(573100,573131):
        x = "HRR" + str(i)
        cmd_string = "gatk MarkDuplicates -I /mnt/bwa-0.7.17/output/"+x+".marked.RG.bam -O /mnt/bwa-0.7.17/output/"+x+".nodup.bam -M /mnt/bwa-0.7.17/output/"+x+".metrics.txt --REMOVE_DUPLICATES true --ASSUME_SORTED true"
        print('x：{}'.format(cmd_string))
        print(os.popen(cmd_string).read())
                
gatk()

生成索引文件

samtools index /mnt/bwa-0.7.17/output/HRR573099.nodup.bam
大概会生成HRR573099.nodup.bam.bai
批量：

def gatk():
    import os
    for i in range(573100,573131):
        x = "HRR" + str(i)
        cmd_string = "samtools index /mnt/bwa-0.7.17/output/"+x+".nodup.bam"
        print('x：{}'.format(cmd_string))
        print(os.popen(cmd_string).read())
                
gatk()

BQSR（重新校准基础质量得分）

首先，我们需要先为hg38.fa生成dict文件
gatk CreateSequenceDictionary -R /mnt/bwa-0.7.17/output/samtools-1.9/hg38.fa -O /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/hg38.dict && echo "** dict done **"

BQSR分两步进行，分别使用BaseRelibrator与ApplyBQSR两个工具。

第一步：BaseRelibrator，此工具根据物种的已知sites，生成一个校正质量值所需要的表格文件。用之前indel区域重比对的命令可以运行：
gatk BaseRecalibrator -I /mnt/bwa-0.7.17/output/HRR573099.marked.RG.bam -R /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/Homo_sapiens_assembly38.fasta --known-sites /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz --known-sites /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/dbsnp_146.hg38.vcf.gz --known-sites /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz -O /mnt/bwa-0.7.17/output/HRR573099.recal.table
批量：

def gatk():
    import os
    for i in range(573100,573131):
        x = "HRR" + str(i)
        cmd_string = "gatk BaseRecalibrator -I /mnt/bwa-0.7.17/output/"+x+".marked.RG.bam -R /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/Homo_sapiens_assembly38.fasta --known-sites /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz --known-sites /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/dbsnp_146.hg38.vcf.gz --known-sites /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz -O /mnt/bwa-0.7.17/output/"+x+".recal.table"
        print('x：{}'.format(cmd_string))
        print(os.popen(cmd_string).read())
                
gatk()

第二步：ApplyBQSR，此工具执行BQSR的第二步，具体来说，它会根据第一步中BaseRelibrator工具生成的重新校准表来重新校准输入的BAM的质量，并输出重新校准的BAM文件

gatk ApplyBQSR -I /mnt/bwa-0.7.17/output/HRR573099.nodup.bam -R /mnt/bwa-0.7.17/output/samtools-1.9/hg38.fa --bqsr-recal-file /mnt/bwa-0.7.17/output/HRR573099.recal.table -O /mnt/bwa-0.7.17/output/HRR573099.bqsr.bam
成功
批量：

def gatk():
    import os
    for i in range(573100,573131):
        x = "HRR" + str(i)
        cmd_string = "gatk ApplyBQSR -I /mnt/bwa-0.7.17/output/HRR573099.nodup.bam -R /mnt/bwa-0.7.17/output/samtools-1.9/hg38.fa --bqsr-recal-file /mnt/bwa-0.7.17/output/HRR573099.recal.table -O /mnt/bwa-0.7.17/output/HRR573099.bqsr.bam"
        print('x：{}'.format(cmd_string))
        print(os.popen(cmd_string).read())
                
gatk()

下面，可以对生成的BAM文件进行验证，使用ValidateSamFile工具，若显示未找到错误，就可以进行变异检测了。
gatk ValidateSamFile -I /mnt/bwa-0.7.17/output/HRR573099.bqsr.bam
反馈

……
ERROR: Read name A00783:553:HTTYVDSXY:3:2678:27670:32690, Mate not found for paired read
ERROR: Read name A00783:553:HTTYVDSXY:3:2628:1325:12273, Mate not found for paired read
ERROR: Read name A00783:553:HTTYVDSXY:3:1107:19397:14669, Mate not found for paired read
ERROR: Read name A00783:553:HTTYVDSXY:3:1519:3197:21997, Mate not found for paired read
Maximum output of [100] errors reached.
[Fri Aug 11 15:07:40 UTC 2023] picard.sam.ValidateSamFile done. Elapsed time: 17.45 minutes.
Runtime.totalMemory()=2978480128
To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp
Tool returned:
3

如发现错误，则可以使用FixMateInformation工具进行修复

gatk FixMateInformation -I /mnt/bwa-0.7.17/output/HRR573099.bqsr.bam -O /mnt/bwa-0.7.17/output/HRR573099.bqsr.bam

def gatk():
    import os
    for i in range(573100,573131):
        x = "HRR" + str(i)
        cmd_string = "gatk FixMateInformation -I /mnt/bwa-0.7.17/output/"+x+".bqsr.bam -O /mnt/bwa-0.7.17/output/"+x+".bqsr.bam"
        print('x：{}'.format(cmd_string))
        print(os.popen(cmd_string).read())
                
gatk()

变异检测，使用HaplotypeCaller工具，此步极费时间
gatk HaplotypeCaller -R /mnt/bwa-0.7.17/output/samtools-1.9/hg38.fa -I /mnt/bwa-0.7.17/output/HRR573099.bqsr.bam -O /mnt/bwa-0.7.17/output/HRR573099.g.vcf.gz -ERC GVCF
之前报错：
Runtime.totalMemory()=3149398016 htsjdk.samtools.util.RuntimeIOException: /mnt/bwa-0.7.17/output/HRR573099.bqsr.bam has invalid uncompressedLength: -469995472

后来单独运行3099这一个没有报错，产生

现在有的文件有

（那个vcf.2的可能是多的，其他的样本没有）

生成g.vcf.gz文件之后，就可以使用GenotypeGVCFs工具对多个样本执行联合基因分型，目的是将多个样本合成一个文件

--gatk4合并多个vcf文件为一个maf文件命令
gatk CombineGVCFs -R /mnt/bwa-0.7.17/output/samtools-1.9/hg38.fa --variant /mnt/bwa-0.7.17/output/HRR573129.g.vcf.gz --variant /mnt/bwa-0.7.17/output/HRR573146.g.vcf.gz -O /mnt/bwa-0.7.17/mafOutput/output.maf
先前代码
gatk GenotypeGVCFs -R /mnt/bwa-0.7.17/output/samtools-1.9/hg38.fa -V /mnt/bwa0.7.17/output/HRR573099.output.g.vcf.gz -V /mnt/bwa-0.7.17/output/HRR573208.output.g.vcf.gz -V /mnt/bwa-0.7.17/output/HRR573170.output.g.vcf.gz -V /mnt/bwa6,7.17/output/HRR573185.output.g.vcf,gz -V /mnt/bwa-0,7.17/output/HRR573142.output.g.vcf.gz -V /mnt/bwa-0.7.17/output/HRR573126.output.g.vcf.gz -V /mnt/bwa0.7.17/output/HRR573186.output.g.vcf.gz -V /mnt/bwa-0.7.17/output/HRR573209 .output.g.vcf.gz -V /mnt/bwa-0.7.17/output/HRR573171.output.g.vcf.gz -V /mnt/bwa0.7.17/output/HRR573143.output.g.vcf.gz -V /mnt/bwa-0.7.17/output/HRR573127.output.g.vcf.gz -O /mnt/bwa-0.7.17/output/com.g.vcf.gz
目前可以运行的代码（单个文件）
gatk GenotypeGVCFs -R /mnt/bwa-0.7.17/output/samtools-1.9/hg38.fa -variant /mnt/bwa-0.7.17/output/HRR573099.g.vcf.gz -O /mnt/bwa-0.7.17/output/com.g.vcf.gz

pVACseq
A cancer immunotherapy pipeline for identifying and prioritizing neoantigens from a list of tumor mutations.
一种癌症免疫治疗用于从肿瘤突变列表中识别和优先选择新抗原的方法。
pVACseq是一种结合肿瘤突变和表达数据(DNA-和RNA-Seq)，通过癌症测序(cancer Sequencing, pVACseq)识别个性化变异抗原的癌症免疫治疗管道。它通过使用大量的平行序列数据来预测肿瘤特异性突变多肽(新抗原)，从而使癌症免疫治疗研究成为可能。它被用于检查点治疗反应的研究，以及为个性化癌症疫苗和过继T细胞治疗确定靶点。

教程“VCF vs MAF |变异注释及整理为Maf格式”
GATK Funcotator注释VCF、并自动输出Maf格式的参考代码

提取SNV，使用SelectVariants工具

gatk SelectVariants -R /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/Homo_sapiens_assembly38.fasta -V /mnt/bwa-0.7.17/output/HRR573146.g.vcf.gz --select-type-to-include SNP -O /mnt/bwa-0.7.17/output/HRR573146.snv.g.vcf

提取INDEL，仍然使用SelectVariants工具，只是参数选择INDEL

突变位点质量值重矫正（VQSR）

这一步是对突变位点进行质控，原理是通过比对与已知变异位点的交集，评估各个突变位点的可信度。已知变异集包括：hapmap，omni，1000G，dbsnp。

第一步：SNP的VQSR的过滤
gatk VariantRecalibrator -R /mnt/bwa-0.7.17/output/samtools-1.9/hg38.fa -V /mnt/bwa-0.7.17/output/com.g.vcf.gz --resource hapmap,known=false,training=true,truth=true,prior=15.0:/mnt/bwa-0.7.17/output/hapmap_3.3_grch38_pop_stratified_af.vcf.gz --resource omni,known=false,training=true,truth=false,prior=12.0:/mnt/bwa-0.7.17/output/1000G_omni2.5.hg38.vcf.gz --resource 1000G,known=false,training=true,truth=false,prior=10.0:/mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz --resource dbsnp,known=true,training=false,truth=false,prior=2.0:/mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/dbsnp_146.hg38.vcf.gz -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR -mode SNP -O /mnt/bwa-0.7.17/output/HRR573146.snp.recal --tranches-file /mnt/bwa-0.7.17/output/HRR573146.snp.tranches --rscript-file /mnt/bwa-0.7.17/output/HRR573146.snp.plots.R

报错如下是没有放对应的tbi文件

报错如下是在-an参数所添加的某一个注释信息在vcf文件里面没有，这时候有两种选择，一种是用variant-annotator这个工具在vcf文件中加上这个注释，另一种是把那个-an参数去掉（后来用了GVCF后的com.g.vcf.gz，需要加上-an参数才能用）

INDEL的VQSR的过滤
gatk VariantRecalibrator -R /mnt/bwa-0.7.17/output/samtools-1.9/hg38.fa -V /mnt/bwa-0.7.17/output/com.g.vcf.gz --max-gaussians 4 --resource mills,known=false,training=true,truth=true,prior=12.0:/mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz --resource dbsnp,known=true,training=false,truth=false,prior=2.0:/mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/dbsnp_146.hg38.vcf.gz -an QD -an DP -an MQRankSum -an ReadPosRankSum -an FS -an SOR -mode INDEL -O /mnt/bwa-0.7.17/output/HRR573146.indel.recal --tranches-file /mnt/bwa-0.7.17/output/HRR573146.indel.tranches --rscript-file /mnt/bwa-0.7.17/output/HRR573146.indel.plots.R

ApplyVQSR to SNP

输入文件是由GenotypeGVCFs生成的vcf.gz
gatk ApplyVQSR -R /mnt/bwa-0.7.17/output/samtools-1.9/hg38.fa -V /mnt/bwa-0.7.17/output/com.g.vcf.gz -O /mnt/bwa-0.7.17/output/com.vqsr.snp.vcf.gz --truth-sensitivity-filter-level 99.0 --tranches-file /mnt/bwa-0.7.17/output/HRR573146.snp.tranches --recal-file /mnt/bwa-0.7.17/output/HRR573146.snp.recal -mode SNP
INDEL也是如此
gatk ApplyVQSR -R /mnt/bwa-0.7.17/output/samtools-1.9/hg38.fa -V /mnt/bwa-0.7.17/output/com.g.vcf.gz -O /mnt/bwa-0.7.17/output/com.vqsr.indel.vcf.gz --truth-sensitivity-filter-level 99.0 --tranches-file /mnt/bwa-0.7.17/output/HRR573146.indel.tranches --recal-file /mnt/bwa-0.7.17/output/HRR573146.indel.recal -mode INDEL

SNP文件注释

gatk CNNScoreVariants -V /mnt/bwa-0.7.17/output/com.vqsr.snp.vcf.gz -R /mnt/bwa-0.7.17/output/samtools-1.9/hg38.fa -O /mnt/bwa-0.7.17/output/com.anno.snp.vcf --disable-avx-check true

INDEL文件注释

gatk CNNScoreVariants -V /mnt/bwa-0.7.17/output/HRR573146.indel.g.vcf -R /mnt/bwa-0.7.17/output/samtools-1.9/hg38.fa -O /mnt/bwa-0.7.17/output/HRR573146.anno.indel.vcf

我觉得CNNScoreVariants可能就是在注释。也许不需要annovar，但是CNNScoreVariants一直卡着不动
应该是注释后的文件给maftools
闲鱼上的人说“这个CNNScoreVariants+FilterVariantTranches是适合单样本做打分过滤的，你昨天做的VQSR相当于已经做过打分过滤了，已经可以进行注释了”

使用Annovar

下载Annovar，tar -zxvf annovar.latest.tar.gz

下载数据库

切换到annovar文件夹，输入perl annotate_variation.pl -downdb -webfrom annovar --buildver hg19 refGene humandb，perl annotate_variation.pl -downdb -webfrom annovar --buildver hg19 clinvar_20221231 humandb

vcf格式转换

路径/mnt/bwa-0.7.17/output/annovar

perl convert2annovar.pl -format vcf4 /mnt/bwa-0.7.17/output/com.vqsr.snp.vcf.gz >com.annovar
也不知道这个结果对不对

简单的数据库注释

perl table_annovar.pl HRR573196.annovar humandb --outfile 196 --buildver hg19 --protocol refGene,clinvar_20221231 --operation g,f --vcfinput

右边是annovar的结果，左边是运行日志，貌似没报错

我们也可以利用ANNOVAR的核心程序 annotate_variation.pl，快速简便的完成单一类型的注释
基于基因
annotate_variation.pl -geneanno -dbtype refGene -buildver hg19 example/ex1.avinput humandb/
基于区域
annotate_variation.pl -regionanno -dbtype cytoBand -buildver hg19 example/ex1.avinput humandb/
基于筛选
annotate_variation.pl -filter -dbtype exac03 -buildver hg19 example/ex1.avinput humandb/

查看注释后文件前三行head -n 3 final.hg19_multianno.txt

取前10列，加上Sample这一列，保存为merged_multianno.2.txt.final

library(maftools)
setwd("E:\\外显子\\figure")
var.annovar.maf = annovarToMaf(annovar = "31sample.txt", 
                               Center = 'NA', 
                               refBuild = 'hg19', 
                               tsbCol = 'Sample', 
                               table = 'refGene',
                               sep = "\t")
write.table(var.annovar.maf,file="var_annovar.maf",quote= F,sep="\t",row.names=F)
var_maf = read.maf(maf ="var_annovar.maf",isTCGA = F)
#summary
plotmafSummary(maf = var_maf, rmOutlier = TRUE, addStat = 'median')

多运行几个样本：
gatk GenotypeGVCFs -R /mnt/bwa-0.7.17/output/samtools-1.9/hg38.fa -variant /mnt/bwa-0.7.17/output/HRR573155.g.vcf.gz -O /mnt/bwa-0.7.17/GVCF/HRR573155.g.vcf.gz

注释完的txt修改并拼接，在E:\外显子的“修改.r”

setwd("E:\\外显子\\merge")
list.files(pattern = "\\.txt$")
lf <-list.files(pattern = ".txt$") 
files <- gsub("\\.txt", "", lf)

for (i in seq_along(files)){
  assign(files[i],data<-read.table(lf[i], header=T,sep="\t"))
  data = data[,1:10]
  paths <- lf[i]
  write.table(data,file=paths, sep = "\t", col.names = TRUE,row.names = FALSE,quote = FALSE)
}

for (i in seq_along(files)){
  assign(files[i],data<-read.table(lf[i], header=T,sep="\t"))
  newname = substr(lf[i],1,3)
  data$Sample<-as.numeric(newname)
  paths <- lf[i]
  write.table(data,file=paths, sep = "\t", col.names = TRUE,row.names = FALSE,quote = FALSE)
}

#拼接
for(i in seq_along(files)){
  assign(files[i],data<-read.table(lf[i], header=T,sep="\t"))
  if(i == 1) {out <- data}
  else { out <- rbind(out,data)}
}
write.table(out,file = "31sample.txt", sep = "\t", col.names = TRUE,row.names = FALSE,quote = FALSE)
#test
data<-read.table("155.hg19_multianno.txt", header=T,sep="\t")

由于没有过滤频率、剔除同义突变，保留有害突变，瀑布图里很多黑色，是因为很多突变，下面是姐姐发的代码，我需要下载rmsk数据库，再做一下table_annovar.pl（不用做convert2annovar了）。“1000g2015aug_all,gnomad_genome,exac03”三个数据库是过滤频率的，但是耗时久，所以先rmsk能少很多突变
下载了rmsk后要安装
perl annotate_variation.pl -downdb -buildver hg19 rmsk humandb/

perl convert2annovar.pl -format vcf4 -allsample -withfreq -includeinfo HRR.filter.vcf > HRR.vcf.avinput

perl table_annovar.pl HRR.vcf.avinput 
 /mnt/bwa-0.7.17/output/annovar/humandb/hg38 -buildver hg38 -remove 
 -out 31sample
 -protocol refGene,esp6500siv2_all,1000g2015aug_all,avsnp150,dbnsfp41a,clinvar_20221231,gnomad_genome,exac03,dbscsnv11,cosmic70,rmsk,ensGene,knownGene   -operation  g,f,f,f,f,f,f,f,f,f,r,g,g   -nastring . --otherinfo

refGene的结果里面可以过滤同义突变
1000g2015aug_all,gnomad_genome,exac03,过滤频率
rmsk过滤重复区的变异，注释到Name=XXX的突变都删掉

下载了1000g2015aug_all数据库

def gatk():
    import os
    for i in range(100,145):
        x = "HRR573" + str(i)
        cmd_string = "perl table_annovar.pl "+x+".annovar /mnt/bwa-0.7.17/output/annovar/humandb -buildver hg19 -remove -out "+x+" -protocol refGene,clinvar_20221231,1000g2015aug_all,rmsk   -operation  g,f,f,r -nastring . --otherinfo"
        print('x：{}'.format(cmd_string))
        print(os.popen(cmd_string).read())
                
gatk()

再awk -F "\t" -v OFS="\t" '{ if ($16<0.05 || $16=="."|| $16=="1000g2015aug_all") print }' 102.txt > 102.awk.txt

过滤的不多，再下其他两个数据库
perl annotate_variation.pl -downdb -buildver hg19 -webfrom annovar gnomad_genome humandb/ &

perl annotate_variation.pl -downdb -buildver hg19 -webfrom annovar exac03 humandb/ &

可以把gnomad30_genome的30去掉

def gatk():
    import os
    for i in range(100,145):
        x = "HRR573" + str(i)
        cmd_string = "perl table_annovar.pl "+x+".annovar /mnt/bwa-0.7.17/output/annovar/humandb -buildver hg19 -remove -out "+x+" -protocol refGene,clinvar_20221231,gnomad_genome,1000g2015aug_all,exac03,rmsk   -operation  g,f,f,f,f,r -nastring . --otherinfo"
        cmd = "grep -v \"stream\" "+x+ ".hg19_multianno.txt|grep -v \"UTR\"|grep -v \"ncRNA\"|grep -v \"        synonymous\"|grep -v \"intronic\"|grep -v \"intergenic\"|grep -v \"nonframeshift\"|grep -v \"Name=\" > "+x+".txt"
        print(f'x：{cmd}')
        print(os.popen(cmd).read())
                
gatk()

def gatk():
    import os
    for i in range(100,145):
        x = "HRR573" + str(i)
        cmd_string = "perl table_annovar.pl "+x+".annovar /mnt/bwa-0.7.17/output/annovar/humandb -buildver hg19 -remove -out "+x+" -protocol refGene,clinvar_20221231,gnomad_genome,1000g2015aug_all,exac03,rmsk   -operation  g,f,f,f,f,r -nastring . --otherinfo"
        cmd = "grep -v \"stream\" "+x+ ".hg19_multianno.txt|grep -v \"UTR\"|grep -v \"ncRNA\"|grep -v \"        synonymous\"|grep -v \"intronic\"|grep -v \"intergenic\"|grep -v \"nonframeshift\"|grep -v \"Name=\" > "+x+".txt"
        cmd1 = "awk -F \" \\t\" -v OFS=\"\\t\" \'{ if ($16<0.05 || $16==\".\"|| $16==\"1000g2015aug_all\") print }\'  "+x+".txt >  "+x+".awk.txt"
        print(f'x：{cmd1}')
        print(os.popen(cmd1).read())
                
gatk()

过滤方法还是不好，这是姐姐对我数据的处理方式：

加一个COSMIC数据库
B是保留编码区的突变，并且过滤掉同义突变，C是在B的基础上对人群频率进行过滤，具体是：保留在COSMIC存在的突变，且在1000G频率小于10%，ExAC数据库中频率小于10%;保留在COSMIC和dbSNP中不存在的突变，保留ExAC数据库中频率小于0.1%，在1000G频率小于0.3%；

瀑布图里多了很多突变，但大部分还是黑色
这是她处理的命令：
过滤

input_dir=./anno
output_dir=./output
mkdir ${output_dir}
cat ${input_dir}/${sample}.hg19_multianno.txt|grep -v "intergenic"|grep -v "stream"|grep -v "UTR"|grep -v "ncRNA"|grep -v "     synonymous"|grep -v "intronic"|grep -v "nonframeshift" |awk -F "\t" -v OFS="\t" '{ if (($16<0.005 || $16==".")&&($24<0.005 || $24==".")&&($25<0.005 || $25==".")) print }'|grep -v "Name="|grep -v "enign" > ${output_dir}/${sample}.anno.txt

绘图

cd output
#单独运行
head -1 
rm -f ../all_multianno.txt
for i in *.txt
do 
 sample=`echo $i|awk -F '.' '{print $1}'`
 sed "s/$/\t${sample}/" $i >> ../all_multianno.txt 
done
cd ..
cut -f 1-10,37 all_multianno.txt > all_multianno1.txt
cat head.txt all_multianno.txt > merged_multianno.txt.final
rm -f all_multianno*

2.plot.sh (END)

需要找参考文献看过滤条件，每个疾病的背景都不一样，需要找到研究的那个疾病的相关参考
然后我发现我做‘体细胞突变检测’这步错了，用HaplotypeCaller是检测生殖系的。要明确想要检测的突变类型，一般肿瘤是体细胞突变。文献是用Mutect2，体细胞变异会少很多
以下是文献的处理步骤

我做到ApplyBQSR这一步还是对的，下一步应该用Mutect2而不是gatk
路径/mnt/bwa-0.7.17/output/annovar，gatk.py

def gatk():
    import os
    for i in range(573098,573120):
        x = "HRR" + str(i)
        cmd_string = "gatk Mutect2 -R /mnt/bwa-0.7.17/output/samtools-1.9/hg38.fa --dbsnp /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/dbsnp_146.hg38.vcf --cosmic CosmicCodingMuts.vcf -I /mnt/bwa-0.7.17/output/"+x+".bqsr.bam -O /mnt/bwa-0.7.17/Mutect2/"+x+".vcf.gz --tumor-sample "+x
        print('x：{}'.format(cmd_string))
        print(os.popen(cmd_string).read())
                
gatk()

小插曲：samtools突然没有了，文件夹下所有东西没有了，用conda install -c bioconda samtools下载但路径文件夹里什么也没有（1.5版本而且也没有新文件夹产生），下载了tar.bz2上传到文件夹，移动了hg38.fa和hg38.dict到文件夹下

新mutect2命令

路径/mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0
gatk Mutect2 -R /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/Homo_sapiens_assembly38.fasta -I /mnt/bwa-0.7.17/output/HRR573108.nodup.bam -I /mnt/bwa-0.7.17/output/HRR573109.nodup.bam -tumor-sample HRR573108 -normal-sample HRR573109 -O 0809_somatic.vcf.gz -bamout 0809_tumor_normal.bam（这个结果是不对的）
以下文件没有，不填
-pon resources/chr17_m2pon.vcf.gz
-germline-resource resources/chr17_af-only-gnomad_grch38.vcf.gz
-af-of-alleles-not-in-resource 0.0000025
--disable-read-filter MateOnSameContigOrNoMappedMateReadFilter
-L chr17plus.interval_list

/mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/0809_somatic.vcf.gz

samtools view -h HRR573150.bqsr.bam | more
samtools view xx.bam chr1:10000000| more

认识了bam格式，新的mutect2需要肿瘤和正常的一对输入（所以之前的call不对，重新call），发现都是P和C来自同一个人，认为c是肿瘤p是正常，输入文件nodup.bam，报错“htsjdk.samtools.util.RuntimeIOException: /mnt/bwa-0.7.17/output/HRR573108.nodup.bam has invalid uncompressedLength: -1961771837”，认为之前的bam可能有问题，因为质量很低，使用speedseq工具来纠正bam，然后服务器坏了

150M前的数字的意思：

看到每个病人都有P和C两种数据

去GSA看对P和C的解释，没有找到

可以用speedseq这个整合工具，输入Reads，得到符合要求的bam。或者输入现有bam，跑realign，得到符合要求的bam

sppeedseq命令：在路径/mnt/bwa-0.7.17/output/speedseq-0.1.2/bin输入（需要conda activate python27）

./speedseq realign -t 40 /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/Homo_sapiens_assembly38.fasta /mnt/bwa-0.7.17/output/HRR573109.nodup.bam &

当编译不成功可尝试的解决办法有：
sudo make install
make和make install分开运行
有时候即使报错但是也安装完了，make check试试，然后‘工具名 -h’看有没有帮助文档
以下是配置speedseq.config的情况

按按照INSTALLATION的步骤来，但是一直报错:

这是installation:

make && make install分开运行，可以make但不能make install

Error: pysam is not installed for /usr/bin/python2.7

在bin下找到speedseq，把黄框这几行注释掉

因为conda 环境下无法识别 -c 运行指令，直接注释掉就行了，只要保证那些包安装好，就可以跳过包检查环节。
把speedseq装好了
但是bam结果很小

尝试用bqsr.bam试试

speedseq/bin下面这个bam大小合理，之前那个大小很小的bam是我点开下面两个文件夹看到的

用它们做mutect2
/mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0
gatk Mutect2 -R /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/hg38/Homo_sapiens_assembly38.fasta -I /mnt/bwa-0.7.17/output/speedseq-0.1.2/bin/HRR573108.nodup.realign.bam -I /mnt/bwa-0.7.17/output/speedseq-0.1.2/bin/HRR573109.nodup.realign.bam -tumor-sample HRR573108 -normal-sample HRR573109 -O 0809_somatic.vcf.gz -bamout 0809_tumor_normal.bam -native-pair-hmm-threads 40 &

gatk4不能多线程

运行完mutect2，看看多少行：
gunzip -c 0809_somatic.vcf.gz | wc -l

试试NeoPredPipe

/mnt/bwa-0.7.17/output/NeoPredPipe
python NeoPredPipe.py -I /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/vcf -H /mnt/bwa-0.7.17/output/arcasHLA/arcasHLA/HLA/hlatypes.txt -o ./ -n TestRun -c 1 2 -E 8 9 10

之前用的hla是scRNA-seq来的
报‘No input vcf files detected. Perhaps they are compressed’，查看了0809_somatic.vcf是有一万多行的

示例数据python NeoPredPipe.py -I ./Example/input_vcfs -H ./Example/HLAtypes/hlatypes.txt -o ./ -n TestRun -c 1 2 -E 8 9 10

我把vcf和hla都改成了和示例一样的格式，但是报错和之前一样，找不到vcf文件，vcf也不是压缩的。预测的hla在/mnt/bwa-0.7.17/output/arcasHLA/arcasHLA/HLA/*genotype.json

python NeoPredPipe.py -I /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/0809_2_somatic.vcf -H /mnt/bwa-0.7.17/output/arcasHLA/arcasHLA/HLA/hlatypes.txt -o ./ -n TestRun -c 1 2 -E 8 9 1

这段代码，应是在指定路径中搜索vcf后缀的文件

python NeoPredPipe.py -I /mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0/vcf -H /mnt/bwa-0.7.17/output/arcasHLA/arcasHLA/HLA/hlatypes.txt -o ./ -n TestRun -c 1 2 -E 8 9 1

放了8个样本，但是table是空的

好像是netMHCpan的问题
把vcf里文件的名字改成和hla里txt的一样

查看hla确实是24:555

然后我把hla已经改成和示例一样的形式（hla_a_11_303），一样。尝试了各种形式的hla输入格式，HLA-A24:555等都不行。有没有可能arcasHLA得到的有些hla，根本是不存在的？但我发现库里有的hla也识别不了（库：https://github.com/YanqiangLi/NeoPredPipe/blob/master/netMHCpanAlleles.txt），去找arcasHLA调用的库，是IMGT/HLA database。如果用示例里的hla，不报错，但是结果是空的。然后把所有hla换成库里有的，搜不到的就NA，出的报错是示例数据出的报错：

我怀疑这个工具挂了，github上3年前的bug都还没修复。
然后我猜想，是不是vcf文件不合适，然后我看了运行的log，写的是invalid record in VCF file: the GT specifier is not present in the FORMAT string。vcf文件中GT代表基因型，是不是mutect2加一些参数，可能会出来GT的结果。
然后老师说列中圈的就是GT

老师说记得这个程序有设置选哪些列的参数

正确命令

python NeoPredPipe.py -I /mnt_wangnan/bwa/wes/gatk4/gatk-4.0.6.0/vcf -H /mnt_wangnan/bwa/output/arcasHLA/arcasHLA/HLA/hlatypes.txt -o ./ -n TestRun -c 0 -E 8 9 10

发现netMHCpan-4.0是可以用的，4.1是最新版还没经过测试

不能在根目录下新建文件夹

但结果还是0

去看运行日志

filefasta结果正常

又换了netMHCpan版本

运行了8个小时第三个样本还没运行完
但有一些中间文件

外显子测序分析