最近很多人在问我,做生物信息分析,具体是在做哪方面的工作呢?难不难?每天的工作内容是什么?容易上手吗?
先来告诉大家生物信息分析到底难不难。就像我问大家做生物实验难不难,对于做了很久实验的你们来说肯定要告诉我不难,同样的我也要告诉大家生物信息分析也不难;但就像做实验我们需要投入时间精力一样,做生物信息分析或者想要把生物信息分析做好,我们也需要投入大量的时间精力。
那我们做生物信息分析,每天都在做什么呢?
我们先来了解一下生物信息学。很多人认为生物信息学是一门学科,当然是没有错的,但是这个学科里也是有很多分支的。我们有人在做人相关的研究,这个也是大量生物信息学者从事的领域;但除了人有关的研究外,有研究者在做水稻的研究,也有研究者在做大鼠小鼠等各种模式生物的研究。拿人的研究来说,有人通过芯片来做研究,也有人通过测序来做研究。无论是芯片也好,测序也好,都是一种手段,最终都是为了解决生物学家所面对的生物学问题。我们再拿测序来说,有人做DNA-seq,有人做RNA-seq,有人做ChIP-seq,甚至最近已经有了蛋白质-seq。现在生物信息学工作基本也是由着几部分组成,DNA-seq的数据分析处理,RNA-seq的数据分析处理以及ChIP-seq的数据分析处理。
DNA-seq我们在做什么?以及意义是什么?
做DNA-seq最终我们能得到什么呢?拿肿瘤来说,我们可以找到肿瘤中特异发生的SNV(single nucleotide variant),indel(insertion和deletion),SV(结构变异),CNV(copy number variant)等,通过DNA水平的突变分析,我们可以找出哪些突变会对这类肿瘤的发生发展起到了重要的作用,从而了解这类肿瘤的发病机理,以及后续对此类肿瘤有针对性的治疗以及预后提供指导意义。这个指导意义如何体现呢?比如针对某一种癌症现在已经研制出了几种药物,医生如何决定给病人用哪种药物呢?最初始的时候医生只能让病人一个药物一个药物的去尝试,我们都知道癌症的药物都非常贵,一种药物都不一定是一个家庭可以负担的起的,何况把多个药物都尝试一遍呢,不仅花费昂贵,而且也很耽误病人宝贵的治疗时间。但是我们通过DNA-seq发现发生不同突变对不同的药物有较好或差的响应,这样我们就可以通过DNA-seq发现的DNA水平的改变来决定使用或者不适用此种药物。
RNA-seq我们在做什么呢?以及意义是什么?
RNA-seq比较常见的分析方法是拿一批肿瘤样本和一批和肿瘤邻近的正常样本做RNA-seq,通过RNA-seq的data analysis,来确定肿瘤样本和正常样本中表达显著差异的基因,再接着对这些基因做进一步的分析,从而找出在表达水平上哪些基因在对这个疾病起作用,以及是如何起作用的。
当然这个是最基本的RNA-seq分析的流程思路,除此之外,还有lncRNA,circRNA等的分析。
ChIP-seq在做什么?以及意义是什么?
我们知道DNA-seq测的是DNA,RNA-seq测的是RNA,那ChIP-seq测的是什么呢?ChIP,是Chromatin Immunoprecipitation,染色质免疫共沉淀。通过染色质免疫共沉淀的方法特异性的富集与目的蛋白结合的DNA片度,接着再对与目的蛋白结合的DNA进行测序分析,所以我们拿到的数据是DNA的数据。主要用于转录因子、组蛋白等的分析。
这是从几个大的方面对生物信息学的工作做了介绍,那具体是怎么操作的呢?我们拿DNA-seq做个例子介绍一下。
现在市面上的测序仪很多,每个都有各自的特点,但每个测序仪下机后我们拿到的数据都是fastq格式的文件,格式如下所示:
@NCCL-08c6aae1-4fde-46d2-822e-e527b47fe184/1
GGATTACTTGAGCCCAGGAGTTTGAGGTTGTAGTGAGCTGTGTATNNTTGTNTCAAGTAAGAATTTNTNAGTTTTTATTATAANANAATTAGCACAATNN
+
ccccchhhhhhhhhhhhhhghhgggghfghhhhggghhhhhhghgBB[[cdB[[fghhhhhhhghgB[B[[cffgghhhhhghB[B[[dghhhhhhhdBB
##其中每条reads分四行,第一行是意@开头的序列ID信息;第二行是序列的碱基;第三行由“+”开始;第四行为碱基对应的质量,和第二行是一一对应关系。
了解了最初始的文件,我们来看看接下来做什么分析
1.QC
我们拿到了原始数据,需要对我们的数据做质量控制,也就是QC,其中比较常用的简单软件fastqc,fastqc把结果直接做成了可视化的图,其中每个base的质量如下图所示:
2.比对,sort,dedup,realignment,BQSR
我们需要把我们的数据比对到参考基因组上,所以我们需要先把参考基因组信息下载下来,如果我们测的是人的数据,虽然现在已经有了hg38版本的参考基因组,但比较常用的还是hg19,所以我们在选择参考基因组的版本的时候也要注意。目前比较常用的做DNA比对的软件是bwa,用法如下:
bwa mem -t 4 -M -R "\@RG\tID:{library}\tLB:{library}\tPL:Illumina\tPU:{sample}\tSM:{sample}\" ref.fa read1.fastq read2.fastq > test.sam 2> ./test-pe.log
做完比对,我们需要把输出的sam文件转化为bam文件,因为sam文件太大,bam文件相对大小较小;转化为bam文件后需要对bam文件做排序;再接下来需要对排序后的bam文件做去重。用samtools:
samtools view -hbS test.sam -o test.bam
samtools sort test.bam -o test.sorted.bam ##对bam文件做排序
samtools rmdup -S test.sorted.bam test.sorted.rmdup.bam
samtools index test.sorted.rmdup.bam
接着用GATK对bam文件做重新比对。如下所示:
java -jar GenomeAnalysisTK.jar \
-R hg19.fa \
-T RealignerTargetCreator \
-I test.sorted.rmdup.bam \
-o hg19.intervals \
-known Mills_and_1000G_gold_standard.indels.hg19.vcf \
-known 1000G_phase1.indels.hg19.vcf \
java -jar GenomeAnalysisTK.jar \
-R hg19.fa \
-T IndelRealigner \
-targetIntervals hg19.intervals \
-I test.sorted.rmdup.bam \
-o test.sorted.rmdup.realn.bam \
-known Mills_and_1000G_gold_standard.indels.hg19.vcf \
-known 1000G_phase1.indels.hg19.vcf
再接下来需要做BQSR,也就是Base quality score recalibration,步骤如下:
java -jar GenomeAnalysisTK.jar \
-T BaseRecalibrator \
-R hg19.fa \
-I test.sorted.rmdup.realn.bam \
-knownSites dbsnp_137.hg19.vcf \
-knownSites Mills_and_1000G_gold_standard.indels.hg19.vcf \
-knownSites 1000G_phase1.indels.hg19.vcf \
-o test.recal.grp
java -jar GenomeAnalysisTK.jar \
-T PrintReads \
-R hg19.fa \
-I test.sorted.rmdup.realn.bam \
-BQSR test.recal.grp \
-o ChrALL.100.sam.recal.08-3.grp.bam
3.call SNV and INDEL
变异的检测现在也是有比较多的软件,我比较常用的是mutect2,介绍mutect2的用法:
##对于有配对的样本
java -jar GenomeAnalysisTK.jar -T MuTect2 \
-R hg19.fa \
-I:tumor tumor.bam \
-I:normal normal.bam \
--fix_misencoded_quality_scores \
--dbsnp dbsnp_137.hg19.vcf \
--cosmic hg19_cosmic_v54_120711.vcf \
-o tumor_vs_normal_mutect2.vcf
##只有tumor的样本
java -jar GenomeAnalysisTK.jar -T Mutect2 \
-R hg19.fa \
-I sample.bam \
-tumor sample_name \
-O single_sample.vcf.gz
最后再对我们的call出来的SNV和INDEL做注释,那我们的SNV和INDEL分析就完成了。