其他分享
首页 > 其他分享> > 生信 cookbook (持续更新)

生信 cookbook (持续更新)

作者:互联网

目录

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 分别进行过滤,然后合并即可
  1. 体细胞突变
    待更新

  1. 线粒体突变
    待更新


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

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

6. 如何查询 vcf 中突变

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

7. 什么是 PED 文件

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

表型列可选的值有:


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>
  1. Hisat2
  2. 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