生信面试实操分析题¶
一句话概述：面试实操题考察你能否从原始数据出发，独立完成一个完整的分析流程并给出生物学解读，关键是思路清晰、流程规范、结果可靠。
核心知识点速查表¶
| 概念 | 白话解释 |
|---|---|
| 分析思路 | 拿到数据后第一件事：看数据→定方案→跑流程→验证→解读 |
| 可重现性 | 别人拿你的代码+数据能得到相同结果 |
| 质控先行 | 任何分析都从数据质量检查开始 |
| 多重验证 | 用不同方法/工具验证同一个结论 |
| 生物学解读 | 不只是跑出数字，要说明生物学意义 |
一、实操题1：RNA-seq 差异表达分析¶
题目¶
完整解答¶
# ========== 第1步：数据质控 ==========
# 先看看数据长什么样
zcat sample1_R1.fastq.gz | head -8 # 看前2条reads
ls -lh *.fastq.gz # 看文件大小
# 运行FastQC
mkdir -p qc_raw
fastqc *.fastq.gz -o qc_raw/ -t 8 # 原始数据质控报告
# 检查关键指标：
# - 碱基质量：Q30以上？
# - 接头污染：有没有adapter序列？
# - GC分布：是否正态？
# - 重复率：是否过高？
# ========== 第2步：过滤和去接头 ==========
for sample in Control1 Control2 Control3 Treat1 Treat2 Treat3; do
fastp \
-i ${sample}_R1.fastq.gz \ # R1输入
-I ${sample}_R2.fastq.gz \ # R2输入
-o clean/${sample}_R1.clean.fq.gz \ # R1输出
-O clean/${sample}_R2.clean.fq.gz \ # R2输出
-q 20 -l 36 \ # Q20质量+最短36bp
--detect_adapter_for_pe \ # 自动检测接头
-w 4 \ # 4线程
-j clean/${sample}.json \ # JSON报告
-h clean/${sample}.html # HTML报告
done
# ========== 第3步：比对到参考基因组 ==========
# 建STAR索引（只需做一次）
STAR --runMode genomeGenerate \
--genomeDir star_index/ \
--genomeFastaFiles genome.fa \ # 参考基因组
--sjdbGTFfile genes.gtf \ # 基因注释
--runThreadN 8
# 比对
for sample in Control1 Control2 Control3 Treat1 Treat2 Treat3; do
STAR --runThreadN 8 \
--genomeDir star_index/ \
--readFilesIn clean/${sample}_R1.clean.fq.gz \
clean/${sample}_R2.clean.fq.gz \
--readFilesCommand zcat \ # 读压缩文件
--outSAMtype BAM SortedByCoordinate \ # 输出排序BAM
--outFileNamePrefix aligned/${sample}_ \
--quantMode GeneCounts # 同时做基因计数
done
# 检查比对率（应该>80%）
for f in aligned/*_Log.final.out; do
echo "=== $(basename $f) ==="
grep "Uniquely mapped" $f # 唯一比对率
done
# ========== 第4步：基因定量 ==========
# 用featureCounts生成计数矩阵
featureCounts \
-T 8 \ # 8线程
-p --countReadPairs \ # PE模式
-a genes.gtf \ # 基因注释
-o counts.txt \ # 输出文件
aligned/*_Aligned.sortedByCoord.out.bam # 所有BAM文件
# 提取计数矩阵（去掉前面的注释行）
cut -f1,7- counts.txt | tail -n +2 > count_matrix.tsv
# ========== 第5步：DESeq2差异分析 ==========
library(DESeq2)
# 读取计数矩阵
counts <- read.delim("count_matrix.tsv", row.names=1)
# 样本信息
coldata <- data.frame(
condition = factor(c(rep("Control", 3), rep("Treatment", 3)),
levels = c("Control", "Treatment")),
row.names = colnames(counts)
)
# 创建DESeq2对象
dds <- DESeqDataSetFromMatrix(
countData = counts,
colData = coldata,
design = ~ condition
)
# 预过滤（去掉低表达基因）
keep <- rowSums(counts(dds)) >= 10 # 总计数>=10
dds <- dds[keep, ]
# 运行差异分析
dds <- DESeq(dds) # 核心步骤
res <- results(dds, alpha = 0.05) # 提取结果
summary(res) # 摘要
# 结果过滤
sig <- res[which(res$padj < 0.05 & abs(res$log2FoldChange) > 1), ]
sig <- sig[order(sig$padj), ] # 按padj排序
write.csv(as.data.frame(sig), "significant_genes.csv")
cat(sprintf("上调: %d, 下调: %d\n",
sum(sig$log2FoldChange > 0),
sum(sig$log2FoldChange < 0)))
# ========== 第6步：可视化 ==========
# PCA图
vsd <- vst(dds, blind = FALSE)
plotPCA(vsd, intgroup = "condition")
# 火山图
library(EnhancedVolcano)
EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'pvalue',
pCutoff = 0.05,
FCcutoff = 1)
# ========== 第7步：功能富集分析 ==========
library(clusterProfiler)
library(org.Hs.eg.db)
# 基因ID转换
gene_list <- rownames(sig)
entrez <- bitr(gene_list, fromType="SYMBOL",
toType="ENTREZID",
OrgDb=org.Hs.eg.db)
# GO富集
go_result <- enrichGO(
gene = entrez$ENTREZID,
OrgDb = org.Hs.eg.db,
ont = "BP", # 生物过程
pAdjustMethod = "BH",
pvalueCutoff = 0.05
)
dotplot(go_result, showCategory=20)
# KEGG富集
kegg_result <- enrichKEGG(
gene = entrez$ENTREZID,
organism = 'hsa',
pvalueCutoff = 0.05
)
dotplot(kegg_result, showCategory=15)
二、实操题2：变异检测与注释¶
题目¶
完整解答¶
# ========== 变异检测流程 ==========
# 第1步：确认BAM文件质量
samtools flagstat dedup.bam # 比对统计
samtools coverage dedup.bam | head # 覆盖度
picard CollectWgsMetrics \
I=dedup.bam R=ref.fa O=wgs_metrics.txt # WGS指标
# 第2步：BQSR（碱基质量校正）
gatk BaseRecalibrator \
-R reference.fa \
-I dedup.bam \
--known-sites dbsnp.vcf.gz \ # 已知变异位点
--known-sites Mills_and_1000G.vcf.gz \
-O recal_data.table
gatk ApplyBQSR \
-R reference.fa \
-I dedup.bam \
--bqsr-recal-file recal_data.table \
-O recal.bam
# 第3步：HaplotypeCaller变异检测
gatk HaplotypeCaller \
-R reference.fa \
-I recal.bam \
-L targets.bed \ # WES目标区域
-O raw_variants.g.vcf.gz \
-ERC GVCF # 输出GVCF
gatk GenotypeGVCFs \
-R reference.fa \
-V raw_variants.g.vcf.gz \
-O genotyped.vcf.gz
# 第4步：变异过滤
# SNP过滤
gatk SelectVariants -V genotyped.vcf.gz --select-type-to-include SNP -O snps.vcf.gz
gatk VariantFiltration -V snps.vcf.gz \
--filter-expression "QD < 2.0" --filter-name "LowQD" \
--filter-expression "FS > 60.0" --filter-name "HighFS" \
--filter-expression "MQ < 40.0" --filter-name "LowMQ" \
-O filtered_snps.vcf.gz
# INDEL过滤
gatk SelectVariants -V genotyped.vcf.gz --select-type-to-include INDEL -O indels.vcf.gz
gatk VariantFiltration -V indels.vcf.gz \
--filter-expression "QD < 2.0" --filter-name "LowQD" \
--filter-expression "FS > 200.0" --filter-name "HighFS" \
-O filtered_indels.vcf.gz
# 合并
gatk MergeVcfs -I filtered_snps.vcf.gz -I filtered_indels.vcf.gz \
-O filtered_variants.vcf.gz
# 第5步：变异注释
# ANNOVAR注释
perl table_annovar.pl filtered_variants.vcf.gz humandb/ \
-buildver hg38 \
-out annotated \
-protocol refGene,gnomad40_genome,clinvar_20240917,dbnsfp42a \
-operation g,f,f,f \
-nastring . -vcfinput
# 第6步：筛选致病变异
import pandas as pd
variants = pd.read_csv("annotated.hg38_multianno.txt", sep="\t")
print(f"总变异数: {len(variants)}")
# 频率过滤
variants['gnomAD_AF'] = pd.to_numeric(
variants['gnomAD_genome_ALL'], errors='coerce').fillna(0)
rare = variants[variants['gnomAD_AF'] < 0.01] # 罕见变异
print(f"罕见变异(AF<1%): {len(rare)}")
# 功能过滤
functional_types = ['nonsynonymous SNV', 'stopgain', 'frameshift deletion',
'frameshift insertion', 'splicing']
functional = rare[rare['ExonicFunc.refGene'].isin(functional_types)]
print(f"功能变异: {len(functional)}")
# ClinVar注释
pathogenic = functional[
functional['CLNSIG'].str.contains('Pathogenic', na=False)
]
print(f"ClinVar致病: {len(pathogenic)}")
# CADD评分过滤
functional['CADD_phred'] = pd.to_numeric(
functional['CADD_phred'], errors='coerce')
high_impact = functional[functional['CADD_phred'] > 20] # CADD>20有害
print(f"CADD>20: {len(high_impact)}")
# 输出候选致病变异
candidates = high_impact[['Chr', 'Start', 'End', 'Ref', 'Alt',
'Gene.refGene', 'ExonicFunc.refGene',
'gnomAD_AF', 'CADD_phred', 'CLNSIG']]
candidates.to_csv("candidate_pathogenic.csv", index=False)
print(f"\n候选致病变异: {len(candidates)}")
print(candidates.head(10))
三、实操题3：宏基因组分析¶
题目¶
完整解答¶
# ========== 宏基因组分析流程 ==========
# 第1步：质控+去宿主
for sample in $(cat sample_list.txt); do
# 质控
fastp -i ${sample}_R1.fq.gz -I ${sample}_R2.fq.gz \
-o clean/${sample}_R1.fq.gz -O clean/${sample}_R2.fq.gz \
-q 20 -l 50 -w 4 --detect_adapter_for_pe
# 去宿主序列
bowtie2 -x human_genome \
-1 clean/${sample}_R1.fq.gz -2 clean/${sample}_R2.fq.gz \
--un-conc-gz host_removed/${sample}_R%.fq.gz \
-S /dev/null --threads 8 # 比对结果丢弃
done
# 第2步：物种注释
for sample in $(cat sample_list.txt); do
metaphlan host_removed/${sample}_R1.fq.gz,host_removed/${sample}_R2.fq.gz \
--input_type fastq \
--nproc 8 \
-o metaphlan/${sample}_profile.txt \
--bowtie2db metaphlan_db/
done
# 合并物种表
merge_metaphlan_tables.py metaphlan/*_profile.txt > merged_abundance.tsv
# 第3步：多样性分析 + 差异分析
library(vegan)
library(ggplot2)
# 读取物种丰度表
abundance <- read.delim("merged_abundance.tsv", row.names=1)
abundance <- abundance[grep("s__", rownames(abundance)), ] # 只保留物种级别
# 转置（行=样本，列=物种）
abund_t <- t(abundance)
# Alpha多样性
shannon <- diversity(abund_t, index="shannon") # Shannon多样性
richness <- specnumber(abund_t) # 物种丰富度
# 分组比较
groups <- factor(c(rep("Disease",5), rep("Control",5)))
wilcox.test(shannon ~ groups) # Wilcoxon检验
# Beta多样性 - PCoA
bc_dist <- vegdist(abund_t, method="bray") # Bray-Curtis距离
pcoa_res <- cmdscale(bc_dist, eig=TRUE, k=2) # PCoA
# PERMANOVA（检验组间差异是否显著）
adonis2(bc_dist ~ groups, permutations=999) # p<0.05=显著
# 差异物种分析（LEfSe或Wilcoxon）
for(sp in colnames(abund_t)) {
p <- wilcox.test(abund_t[,sp] ~ groups)$p.value
if(p < 0.05) {
cat(sprintf("%s: p=%.4f\n", sp, p))
}
}
四、面试实操题回答框架¶
无论什么实操题，都按这个框架回答：
1. 确认数据（30秒）
"首先我会检查数据格式和质量"
- 文件格式是什么？数据量多大？
- 有没有样本信息表？
2. 说出分析思路（1分钟）
"我的分析分为X步"
- 画一个简单的流程图
- 每一步用什么工具，为什么
3. 关键步骤详述（2-3分钟）
"最关键的步骤是XXX"
- 说出核心参数和阈值
- 说出为什么选这个参数
4. 质控和验证（30秒）
"在每一步我会检查"
- 比对率、重复率、覆盖度
- 结果的合理性（与文献对比）
5. 结果解读（1分钟）
"最终我发现了"
- 用数字说话
- 给出生物学解释
- 提出局限性和改进方向
五、常见报错与解决¶
| 问题 | 原因 | 解决方法 |
|---|---|---|
| 比对率低(<50%) | 参考基因组不对或污染 | 用FastQ Screen检查物种来源 |
| 差异基因为0 | 过滤太严或样本分组错误 | 检查PCA确认分组，放宽阈值 |
| GO富集无显著结果 | 基因列表太短或ID转换失败 | 用GSEA代替ORA，检查ID格式 |
| 变异数异常多/少 | 过滤参数不合适 | 对比dbSNP已知变异数做sanity check |
| 内存不足 | 数据太大 | 分染色体处理，增加内存 |
六、面试高频问题¶
Q1：拿到一份数据你第一步做什么？
第一步永远是看数据：(1) 看文件格式和大小；(2) 看样本信息表（有几组，每组几个重复）；(3) 快速质控（FastQC/MultiQC）；(4) 检查数据完整性（文件是否损坏，reads数是否合理）。不看数据直接跑流程是最大的错误。
Q2：你的分析结果和文献不一致怎么办？
(1) 先检查自己的分析是否有误（参数、版本、过滤条件）；(2) 看是不是样本差异（物种、组织、疾病阶段不同）；(3) 用不同方法重新分析看结果是否稳健；(4) 如果确认没问题，不一致本身也是有价值的发现，可以讨论原因。
Q3：如何保证分析的可重现性？
(1) 记录所有软件版本（conda list > environment.txt）；(2) 用Snakemake/Nextflow管理流程；(3) 代码放Git版本控制；(4) 用Docker/Singularity容器化环境；(5) 写清楚README和方法描述。
七、速查表¶
# === 生信实操分析速查 ===
# RNA-seq标准流程
fastp → STAR → featureCounts → DESeq2 → ClusterProfiler
# WGS/WES变异检测
fastp → BWA-MEM → MarkDuplicates → BQSR → HaplotypeCaller → VQSR/硬过滤
# 宏基因组标准流程
fastp → 去宿主(BWA/Bowtie2) → MetaPhlAn(物种) → HUMAnN(功能) → LEfSe(差异)
# 关键质控指标
# 比对率 > 80%（RNA-seq > 70%）
# 重复率 < 30%（RNA-seq可能更高）
# Q30碱基比例 > 85%
# 外显子比例 > 60%（WES）
# 差异分析阈值
# |log2FC| > 1 且 padj < 0.05（RNA-seq）
# gnomAD AF < 0.01（罕见变异）
# LDA > 2（LEfSe）
# 面试回答框架
# 1. 确认数据 → 2. 分析思路 → 3. 关键步骤 → 4. 质控验证 → 5. 结果解读
参考资料：GATK Best Practices 2024、RNA-seq分析标准流程、MetaPhlAn4文档、DESeq2 Vignette