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

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

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

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

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

---

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

• 数据预处理：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

#批量执行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这一个没有报错，产生

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

gatk GenotypeGVCFs -R /home/eric/ncbi/hg38/chroms/hg38.fa -V /mnt/d/wes/caller/output.g.vcf.gz -O /mnt/d/wes/caller/output.vcf.gz