一望蜡、外顯子測序分析可以得到以下結(jié)果：

• 變異位點(diǎn)的檢測和注釋
• 基因突變的檢測和注釋
• 拼接剪接位點(diǎn)的檢測和注釋
• 基因外顯子區(qū)域的覆蓋度分析
• 基因外顯子區(qū)域的差異表達(dá)分析
• 基因功能富集分析

二阁猜、利用外顯子測序數(shù)據(jù)得到腫瘤特異性靶標(biāo)的流程一般包括以下幾個(gè)步驟：

• 數(shù)據(jù)預(yù)處理：包括測序數(shù)據(jù)的質(zhì)量控制蜜猾、去除低質(zhì)量序列、比對到基因組參考序列等柿菩。
• 變異檢測：利用不同的軟件對測序數(shù)據(jù)進(jìn)行變異檢測戚嗅，篩選出與腫瘤相關(guān)的變異位點(diǎn)。
• 變異注釋：對變異位點(diǎn)進(jìn)行注釋，包括變異的類型懦胞、位置替久、影響等信息。
• 功能分析：對注釋的結(jié)果進(jìn)行功能富集分析躏尉，篩選出與腫瘤相關(guān)的基因和通路蚯根。
• 靶標(biāo)篩選：根據(jù)功能分析的結(jié)果，篩選出與腫瘤相關(guān)的靶標(biāo)胀糜，進(jìn)行進(jìn)一步的實(shí)驗(yàn)驗(yàn)證颅拦。

具體流程還需要根據(jù)具體的實(shí)驗(yàn)設(shè)計(jì)和數(shù)據(jù)分析要求進(jìn)行調(diào)整。

外顯子測序數(shù)據(jù)分析中常用的軟件和工具：

• 數(shù)據(jù)預(yù)處理：Trimmomatic教藻、FastQC距帅、BWA、SAMtools等括堤。
• 變異檢測：GATK碌秸、VarScan、Mutect2悄窃、Strelka等讥电。
• 變異注釋：ANNOVAR、VEP轧抗、SNPEff等恩敌。
• 功能分析：DAVID、GOseq横媚、KEGG等纠炮。
• 靶標(biāo)篩選：Cytoscape、STRING分唾、GeneMANIA等。

需要注意的是狮斗，具體軟件的選擇還要根據(jù)實(shí)驗(yàn)設(shè)計(jì)绽乔、數(shù)據(jù)類型和分析需求等因素進(jìn)行綜合考慮。

基礎(chǔ)命令：
bwa mem hg38.fa HRR573094_f1.fq.gz HRR573094_r2.fq.gz > mem-pe.3094.sam

可添加的參數(shù)：

• -t：指定線程數(shù)碳褒≌墼遥可以根據(jù)計(jì)算機(jī)配置來設(shè)置。比如沙峻，有8個(gè)CPU核心可用睦授，可以設(shè)置為"-t 8"，以加快比對速度摔寨。
• -M：生成CIGAR字符串中的M標(biāo)記去枷，表示匹配和誤配。這個(gè)選項(xiàng)會讓BWA生成一個(gè)SAM文件，其中包含了所有的比對結(jié)果删顶。如果只需要最好的比對結(jié)果竖螃，可以去掉這個(gè)選項(xiàng)。
• -R：指定read group信息逗余。這個(gè)信息對于后續(xù)的數(shù)據(jù)分析很重要特咆，建議填寫÷剂唬可以參考以下格式：-R "@RG\tID:group1\tSM:sample1\tPL:illumina\tPU:unit1"
• -B：指定比對算法腻格。BWA有三種比對算法可以選擇，分別是"mem"啥繁、"bwasw"和"aln"菜职。其中，默認(rèn)的是"mem"算法输虱，適用于短讀長于70bp的情況些楣。

綜上所述，可以使用以下命令進(jìn)行比對：
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

使用循環(huán)語句來批量處理數(shù)據(jù)宪睹。下面是一個(gè)示例腳本愁茁，可以批量處理一個(gè)目錄下的所有雙端測序數(shù)據(jù)：

bash
#!/bin/bash
# 批量比對雙端測序數(shù)據(jù)

# 1. 定義變量
REF="hg38.fa"
THREADS=8
DIR="data"  # 存放數(shù)據(jù)的目錄
OUTDIR="output"  # 存放比對結(jié)果的目錄

# 2. 創(chuàng)建輸出目錄
mkdir -p $OUTDIR

# 3. 處理每個(gè)樣本
for file in $DIR/*.fq.gz
do

# 3.1 獲取文件名和路徑
filename=$(basename "$file")
sample="${filename%.*}"
   
 # 3.2 執(zhí)行比對命令
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 轉(zhuǎn)換為BAM格式并排序
samtools view -b -S $OUTDIR/${sample}.sam | samtools sort -o $OUTDIR/${sample}.bam -

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

上述腳本中，我們首先定義了參考基因組文件和線程數(shù)等變量亭病，然后創(chuàng)建了一個(gè)用于存放比對結(jié)果的目錄鹅很。接下來使用for循環(huán)遍歷指定目錄下的所有雙端測序數(shù)據(jù)，依次執(zhí)行比對命令罪帖、轉(zhuǎn)換為BAM格式并排序促煮、建立索引等操作。

可以將上述腳本保存為.sh文件整袁，然后使用終端運(yùn)行即可菠齿。需要安裝BWA和SAMtools等軟件，并將它們的可執(zhí)行文件添加到系統(tǒng)環(huán)境變量中坐昙，才能在終端中直接運(yùn)行命令绳匀。

在Illumina測序中，每個(gè)測序文庫都有一個(gè)唯一的標(biāo)識符炸客，稱為“read group”疾棵。這個(gè)標(biāo)識符可以用來區(qū)分不同的測序文庫，并且可以提供有關(guān)測序數(shù)據(jù)的更多信息痹仙，例如測序平臺是尔、測序日期、測序文庫類型等开仰。在數(shù)據(jù)分析中拟枚，read group信息可以幫助我們更好地理解和解釋數(shù)據(jù)薪铜。

下面是一個(gè)read group信息的示例：

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

其中，\t表示制表符梨州，\n表示換行符痕囱。各個(gè)字段的含義如下：

• ID： 測序文庫的唯一標(biāo)識符。
• SM： 樣本名稱暴匠，即測序文庫對應(yīng)的樣本名稱鞍恢。
• PL： 測序平臺，例如Illumina每窖、Ion Torrent等帮掉。
• PU： 測序文庫的物理標(biāo)識符，例如測序文庫的barcode或lane信息等窒典。

這個(gè)示例中的read group信息表明這個(gè)測序文庫的唯一標(biāo)識符是group1蟆炊，對應(yīng)的樣本名稱是sample1，測序平臺是Illumina瀑志，物理標(biāo)識符是unit1涩搓。

在Cell Ranger中，使用-R參數(shù)可以指定read group信息劈猪。如果您不確定如何填寫read group信息昧甘，可以參考Illumina官方文檔或者咨詢測序服務(wù)商。

shell的命令搞不定战得。充边。自己寫的python腳本挺好用。常侦。

用bwa比對

#批量執(zhí)行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 轉(zhuǎn)bam

conda install -c bioconda samtools

#批量執(zhí)行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的命令終于能用了浇冰，網(wǎng)上好多教程是錯(cuò)的

排序

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()

本來的報(bào)錯(cuò)：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 --沒有用聋亡，后來把能運(yùn)行的腳本復(fù)制然后放這個(gè)的代碼肘习。

標(biāo)記重復(fù)序列-去重

#批量執(zhí)行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()

創(chuàng)建索引

進(jìn)行sambamba的markdup，就會生成.bai文件坡倔，所以我想是不是已經(jīng)有index了

index區(qū)域重比對漂佩，重比對區(qū)域校正

發(fā)現(xiàn)網(wǎng)上的教程也是五花八門各不相同。致讥。所以就一個(gè)個(gè)試哪個(gè)能用

生成fai文件：
samtools faidx Homo_sapiens_assembly38.fasta

太久沒wes了其實(shí)工作路徑在/mnt/bwa-0.7.17/wes/gatk4/gatk-4.0.6.0

indel區(qū)域重比對

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

出現(xiàn)報(bào)錯(cuò)“A USER ERROR has occurred: Number of read groups must be >= 1, but is 0”原因是沒有read group分組信息仅仆。
教程：https://uteric.github.io/WES/wes01/
理應(yīng)在bwa的時(shí)候就添加分組信息去跑（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補(bǔ)救一下，示例：
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

用循環(huán)運(yùn)行：
/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（重新校準(zhǔn)基礎(chǔ)質(zhì)量得分）

首先咳榜，我們需要先為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分兩步進(jìn)行，分別使用BaseRelibrator與ApplyBQSR兩個(gè)工具爽锥。

第一步：BaseRelibrator涌韩，此工具根據(jù)物種的已知sites，生成一個(gè)校正質(zhì)量值所需要的表格文件氯夷。用之前indel區(qū)域重比對的命令可以運(yùn)行：
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臣樱，此工具執(zhí)行BQSR的第二步，具體來說腮考，它會根據(jù)第一步中BaseRelibrator工具生成的重新校準(zhǔn)表來重新校準(zhǔn)輸入的BAM的質(zhì)量雇毫，并輸出重新校準(zhǔn)的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文件進(jìn)行驗(yàn)證踩蔚，使用ValidateSamFile工具棚放，若顯示未找到錯(cuò)誤，就可以進(jìn)行變異檢測了馅闽。
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

如發(fā)現(xiàn)錯(cuò)誤飘蚯，則可以使用FixMateInformation工具進(jìn)行修復(fù)

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工具福也，此步極費(fèi)時(shí)間
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
之前報(bào)錯(cuò)：
Runtime.totalMemory()=3149398016 htsjdk.samtools.util.RuntimeIOException: /mnt/bwa-0.7.17/output/HRR573099.bqsr.bam has invalid uncompressedLength: -469995472

后來單獨(dú)運(yùn)行3099這一個(gè)沒有報(bào)錯(cuò)局骤，產(chǎn)生

現(xiàn)在有的文件有

（那個(gè)vcf.2的可能是多的，其他的樣本沒有）

生成g.vcf.gz文件之后拟杉，就可以使用GenotypeGVCFs工具對多個(gè)樣本執(zhí)行聯(lián)合基因分型庄涡，目的是將多個(gè)樣本合成一個(gè)文件

--gatk4合并多個(gè)vcf文件為一個(gè)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
目前可以運(yùn)行的代碼（單個(gè)文件）
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.
一種癌癥免疫治療用于從腫瘤突變列表中識別和優(yōu)先選擇新抗原的方法。
pVACseq是一種結(jié)合腫瘤突變和表達(dá)數(shù)據(jù)(DNA-和RNA-Seq)搬设，通過癌癥測序(cancer Sequencing, pVACseq)識別個(gè)性化變異抗原的癌癥免疫治療管道穴店。它通過使用大量的平行序列數(shù)據(jù)來預(yù)測腫瘤特異性突變多肽(新抗原)，從而使癌癥免疫治療研究成為可能拿穴。它被用于檢查點(diǎn)治療反應(yīng)的研究泣洞，以及為個(gè)性化癌癥疫苗和過繼T細(xì)胞治療確定靶點(diǎn)。

教程“VCF vs MAF |變異注釋及整理為Maf格式”
GATK Funcotator注釋VCF默色、并自動(dòng)輸出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工具腿宰，只是參數(shù)選擇INDEL

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 INDEL -O /mnt/bwa-0.7.17/output/HRR573146.indel.g.vcf

突變位點(diǎn)質(zhì)量值重矯正（VQSR）

這一步是對突變位點(diǎn)進(jìn)行質(zhì)控呕诉，原理是通過比對與已知變異位點(diǎn)的交集，評估各個(gè)突變位點(diǎn)的可信度吃度。已知變異集包括：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

報(bào)錯(cuò)如下是沒有放對應(yīng)的tbi文件

報(bào)錯(cuò)如下是在-an參數(shù)所添加的某一個(gè)注釋信息在vcf文件里面沒有英遭，這時(shí)候有兩種選擇，一種是用variant-annotator這個(gè)工具在vcf文件中加上這個(gè)注釋亦渗，另一種是把那個(gè)-an參數(shù)去掉（后來用了GVCF后的com.g.vcf.gz挖诸，需要加上-an參數(shù)才能用）

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一直卡著不動(dòng)
應(yīng)該是注釋后的文件給maftools
閑魚上的人說“這個(gè)CNNScoreVariants+FilterVariantTranches是適合單樣本做打分過濾的多律，你昨天做的VQSR相當(dāng)于已經(jīng)做過打分過濾了，已經(jīng)可以進(jìn)行注釋了”

使用Annovar

下載Annovar搂蜓，tar -zxvf annovar.latest.tar.gz

下載數(shù)據(jù)庫

切換到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格式轉(zhuǎn)換

路徑/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
也不知道這個(gè)結(jié)果對不對

簡單的數(shù)據(jù)庫注釋

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

右邊是annovar的結(jié)果富纸，左邊是運(yùn)行日志卢佣，貌似沒報(bào)錯(cuò)

我們也可以利用ANNOVAR的核心程序 annotate_variation.pl挤土，快速簡便的完成單一類型的注釋
基于基因
annotate_variation.pl -geneanno -dbtype refGene -buildver hg19 example/ex1.avinput humandb/
基于區(qū)域
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')

多運(yùn)行幾個(gè)樣本：
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")

由于沒有過濾頻率淆院、剔除同義突變此再，保留有害突變昔搂，瀑布圖里很多黑色，是因?yàn)楹芏嗤蛔兪淠矗旅媸墙憬惆l(fā)的代碼摘符，我需要下載rmsk數(shù)據(jù)庫，再做一下table_annovar.pl（不用做convert2annovar了）策吠」淇悖“1000g2015aug_all,gnomad_genome,exac03”三個(gè)數(shù)據(jù)庫是過濾頻率的，但是耗時(shí)久猴抹，所以先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的結(jié)果里面可以過濾同義突變
1000g2015aug_all,gnomad_genome,exac03,過濾頻率
rmsk過濾重復(fù)區(qū)的變異带族，注釋到Name=XXX的突變都刪掉

下載了1000g2015aug_all數(shù)據(jù)庫

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()

下面的代碼也是在linux上運(yùn)行的，要?jiǎng)h一些突變
grep -v "stream" HRR573143.hg19_multianno.txt|grep -v "UTR"|grep -v "ncRNA"|grep -v " synonymous"|grep -v "intronic"|grep -v "intergenic"|grep -v "nonframeshift"|grep -v "Name=" > out.txt
這個(gè)grep -v " synonymous"里面的空的是tab蟀给，用ctrl+V+i鍵打出來蝙砌，這命令真是來之不易

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

過濾的不多，再下其他兩個(gè)數(shù)據(jù)庫
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()

過濾方法還是不好跋理，這是姐姐對我數(shù)據(jù)的處理方式：

加一個(gè)COSMIC數(shù)據(jù)庫
B是保留編碼區(qū)的突變择克，并且過濾掉同義突變，C是在B的基礎(chǔ)上對人群頻率進(jìn)行過濾前普，具體是：保留在COSMIC存在的突變肚邢，且在1000G頻率小于10%，ExAC數(shù)據(jù)庫中頻率小于10%;保留在COSMIC和dbSNP中不存在的突變汁政，保留ExAC數(shù)據(jù)庫中頻率小于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
#單獨(dú)運(yùn)行
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)

需要找參考文獻(xiàn)看過濾條件勺鸦，每個(gè)疾病的背景都不一樣，需要找到研究的那個(gè)疾病的相關(guān)參考
然后我發(fā)現(xiàn)我做‘體細(xì)胞突變檢測’這步錯(cuò)了目木，用HaplotypeCaller是檢測生殖系的换途。要明確想要檢測的突變類型，一般腫瘤是體細(xì)胞突變刽射。文獻(xiàn)是用Mutect2军拟，體細(xì)胞變異會少很多
以下是文獻(xiàn)的處理步驟

我做到ApplyBQSR這一步還是對的，下一步應(yīng)該用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版本而且也沒有新文件夾產(chǎn)生）摹恰，下載了tar.bz2上傳到文件夾辫继，移動(dòng)了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（這個(gè)結(jié)果是不對的）
以下文件沒有，不填
-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

認(rèn)識了bam格式俗慈，新的mutect2需要腫瘤和正常的一對輸入（所以之前的call不對姑宽，重新call），發(fā)現(xiàn)都是P和C來自同一個(gè)人闺阱，認(rèn)為c是腫瘤p是正常炮车，輸入文件nodup.bam，報(bào)錯(cuò)“htsjdk.samtools.util.RuntimeIOException: /mnt/bwa-0.7.17/output/HRR573108.nodup.bam has invalid uncompressedLength: -1961771837”酣溃，認(rèn)為之前的bam可能有問題瘦穆，因?yàn)橘|(zhì)量很低，使用speedseq工具來糾正bam赊豌，然后服務(wù)器壞了

150M前的數(shù)字的意思：

看到每個(gè)病人都有P和C兩種數(shù)據(jù)

去GSA看對P和C的解釋难审，沒有找到

可以用speedseq這個(gè)整合工具，輸入Reads亿絮，得到符合要求的bam告喊。或者輸入現(xiàn)有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 &

當(dāng)編譯不成功可嘗試的解決辦法有：
sudo make install
make和make install分開運(yùn)行
有時(shí)候即使報(bào)錯(cuò)但是也安裝完了，make check試試蒂萎，然后‘工具名 -h’看有沒有幫助文檔
以下是配置speedseq.config的情況

按按照INSTALLATION的步驟來秆吵，但是一直報(bào)錯(cuò):

這是installation:

make && make install分開運(yùn)行，可以make但不能make install

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

因?yàn)閏onda 環(huán)境下無法識別 -c 運(yùn)行指令纳寂，直接注釋掉就行了主穗，只要保證那些包安裝好，就可以跳過包檢查環(huán)節(jié)毙芜。
把speedseq裝好了
但是bam結(jié)果很小

嘗試用bqsr.bam試試

speedseq/bin下面這個(gè)bam大小合理忽媒，之前那個(gè)大小很小的bam是我點(diǎn)開下面兩個(gè)文件夾看到的

用它們做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不能多線程

運(yùn)行完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來的
報(bào)‘No input vcf files detected. Perhaps they are compressed’腋粥，查看了0809_somatic.vcf是有一萬多行的

示例數(shù)據(jù)python NeoPredPipe.py -I ./Example/input_vcfs -H ./Example/HLAtypes/hlatypes.txt -o ./ -n TestRun -c 1 2 -E 8 9 10

我把vcf和hla都改成了和示例一樣的格式晦雨，但是報(bào)錯(cuò)和之前一樣，找不到vcf文件隘冲，vcf也不是壓縮的闹瞧。預(yù)測的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

這段代碼，應(yīng)是在指定路徑中搜索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個(gè)樣本展辞，但是table是空的

好像是netMHCpan的問題
把vcf里文件的名字改成和hla里txt的一樣

查看hla確實(shí)是24:555

然后我把hla已經(jīng)改成和示例一樣的形式（hla_a_11_303）奥邮，一樣。嘗試了各種形式的hla輸入格式罗珍，HLA-A24:555等都不行漠烧。有沒有可能arcasHLA得到的有些hla，根本是不存在的靡砌？但我發(fā)現(xiàn)庫里有的hla也識別不了（庫：https://github.com/YanqiangLi/NeoPredPipe/blob/master/netMHCpanAlleles.txt）已脓，去找arcasHLA調(diào)用的庫，是IMGT/HLA database通殃。如果用示例里的hla度液，不報(bào)錯(cuò)，但是結(jié)果是空的画舌。然后把所有hla換成庫里有的堕担，搜不到的就NA，出的報(bào)錯(cuò)是示例數(shù)據(jù)出的報(bào)錯(cuò)：

我懷疑這個(gè)工具掛了曲聂，github上3年前的bug都還沒修復(fù)霹购。
然后我猜想，是不是vcf文件不合適朋腋，然后我看了運(yùn)行的log齐疙，寫的是invalid record in VCF file: the GT specifier is not present in the FORMAT string。vcf文件中GT代表基因型旭咽，是不是mutect2加一些參數(shù)贞奋，可能會出來GT的結(jié)果。
然后老師說列中圈的就是GT

老師說記得這個(gè)程序有設(shè)置選哪些列的參數(shù)

正確命令

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

發(fā)現(xiàn)netMHCpan-4.0是可以用的穷绵，4.1是最新版還沒經(jīng)過測試

不能在根目錄下新建文件夾

但結(jié)果還是0

去看運(yùn)行日志

filefasta結(jié)果正常

又換了netMHCpan版本

運(yùn)行了8個(gè)小時(shí)第三個(gè)樣本還沒運(yùn)行完
但有一些中間文件