生信 cookbook (持续更新)

2021-06-01 14:35:49 阅读：242 来源： 互联网

标签：10 STAR vcf -- 更新如何 gatk 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 最佳实践？

胚系突变, 这里为了提升运行速度，省略了 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 分别进行过滤，然后合并即可

体细胞突变
待更新

线粒体突变
待更新

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

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>

Hisat2
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

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

ICode9