生信 cookbook (持续更新)

目录

• 1. 如何对 fastq 数据进行过滤？
• 2. 如何使用 bwa mem 生成 bam
• 3. 如何对 bam 文件统计平均深度、覆盖度、捕获情况等信息
• 4. GATK 最佳实践？
• 5. 如何比较两个 vcf 的差异
• 6. 如何查询 vcf 中突变
• 7. 什么是 PED 文件
• 8. 使用 vep 注释
• 9. 如何进行转录组数据比对
• 10. 如何计算基因表达量

1. 如何对 fastq 数据进行过滤？

推荐使用 fastp

fastp -i <fastq1> -I <fastq2>  -o <fastq1.clean> -O <fastq2.clean> -w <cpus> -j <qc.json> -h <qc.html>

---

2. 如何使用 bwa mem 生成 bam

bwa mem -t <cpus> -M -Y -R "@RG\tID:<name>\tPL:ILLUMINA\tSM:<name>" <ref.fastq> <fastq1> <fastq2>
 | samtools view -@ <cpus> -Sb - > <bam>

---

3. 如何对 bam 文件统计平均深度、覆盖度、捕获情况等信息

捕获测序推荐使用 bamdst，全基因组推荐使用

bamdst -p <probe.bed> -o <outdir> <bam>
mosdepth -t <cpus> -n --by <bed> <out_name> <bam>

---

4. GATK 最佳实践？

1. 胚系突变, 这里为了提升运行速度，省略了 BQSR、VQSR 步骤

gatk MarkDuplicates -I <bam> -O <markdup.bam> -M <metrics.txt>
gatk HaplotypeCaller  -R <ref.fasta> -I <markdup.bam>  -L <chr|bed> -O <g.vcf.gz> -ERC GVCF
gatk GenotypeGVCFs -R <ref.fasta>  -V  <g.vcf.gz> -O <vcf.gz>
gatk  SelectVariants -V <vcf.gz> -select-type SNP -O <snp.vcf>
gatk  SelectVariants -V <vcf.gz> -select-type INDEL -O <indel.vcf>
# 对 snp 和 indel 分别进行过滤，然后合并即可

2. 体细胞突变
   待更新

3. 线粒体突变
   待更新

---

5. 如何比较两个 vcf 的差异

# 需安装 vcftools
vcf-compare <vcf1> <vcf2>

---

6. 如何查询 vcf 中突变

bcftools view -r <chr>:<pos> vcf.gz

---

7. 什么是 PED 文件

ped 文件是对样本之间家系关系的描述，主要有 6 列

• Family ID, 家系编号
• Individual ID, 样本编号
• Paternal ID, 样本父亲编号
• Maternal ID, 样本母亲编号
• Sex (1=male; 2=female; other=unknown)
• Phenotype, 表型

表型列可选的值有：

• -9 missing, 缺失
• 0 missing, 缺失
• 1 unaffected, 未受影响（正常）
• 2 affected, 受影响（突变型）

---

8. 使用 vep 注释

vep -i <vcf> \
            -o <out>
            --cache \
            --species <homo_sapiens> \
            --assembly <eg:GRCH37>  \
            --compress_output gzip \
            --force_overwrite \
            --offline \
            --tab \
            --fork <cpus> \
            --everything \
            --plugin <db>\
            --custom <db>\
            --fasta <ref> --merged --hgvs --no_escape

---

9. 如何进行转录组数据比对

推荐使用 Hisat2, STAR 或者 Tophat2, 人类转录组数据推荐使用 STAR

1. STAR

# 建索引
# 如果是人类基因组，可在STAR 官网下载索引
STAR
    --runThreadN <cpus> \
    --runMode genomeGenerate \
    --genomeDir </path/to/genomeDir> \
    --genomeFastaFiles </path/to/genome/fasta1> \
    --sjdbGTFfile </path/to/annotations.gtf> \
    # --sjdbOverhang ReadLength-1
    
 # 比对 
 STAR
    --runThreadN <cpus>
    --genomeDir </path/to/genomeDir>
    --readFilesIn </path/to/read1> </path/to/read2>
    --outFileNamePrefix <name>

2. Hisat2
3. Tophat2

---

10. 如何计算基因表达量

RPKM, 单端的叫 FPKM:

\[ RPKM = \frac{ExonMappedReads*10^9}{TotalMappedReads/GenomeLength} \]

RPM:

\[ RPM or CPM = \frac{ExonMappedReads*10^6}{TotalMappedReads} \]

TPM:

\[ TPM = \frac{N_i/L_i*10^6}{\Sigma_{i=1}^n \frac{N_i}{L_i}} \]

\(N_i\)为比对到第i个exon的reads数； \(L_i\)为第i个exon的长度

RPKM 可以使用 R 包 DESeq2 计算，TPM 使用 egdeR 计算

---

标签：10,STAR,vcf,--,更新,如何,gatk,cookbook,生信
来源： https://www.cnblogs.com/limon86/p/14837053.html