腫瘤空間異質(zhì)性探究

1、腫瘤空間異質(zhì)性探究腫瘤異質(zhì)性包括空間異質(zhì)性、 時(shí)間異質(zhì)性、 解剖異質(zhì)性、 結(jié)構(gòu)異質(zhì)性、基因異質(zhì)性和功能異質(zhì)性等等腫瘤異質(zhì)性是惡 性腫瘤的特征之一，是指腫瘤在生長(zhǎng)過程中，經(jīng)過多次分裂 增殖，其子細(xì)胞呈現(xiàn)出分子生物學(xué)或基因方面的改變，從而 使腫瘤的生長(zhǎng)速度、侵襲能力、對(duì)藥物的敏感性、預(yù)后等各 方面產(chǎn)生差異。腫瘤異質(zhì)性一直是腫瘤治療的挑戰(zhàn)之一，腫 瘤內(nèi)部不同亞群的細(xì)胞對(duì)藥物敏感性的不同可能會(huì)導(dǎo)致治 療的失敗?，F(xiàn)在主流的探究腫瘤異質(zhì)性的方法是：對(duì)腫瘤病 人的腫瘤組織進(jìn)行不同時(shí)間點(diǎn)取樣對(duì)腫瘤病人的腫瘤組織 不同部分分別取樣取樣后測(cè)序可以是 WES,WGS ，或者靶向 部分特殊基因進(jìn)行高深度測(cè)序。數(shù)據(jù)分析

2、的時(shí)候探究不同測(cè) 序結(jié)果的 somatic mutation 的共有比例，或者部分突變的 allele frequency 的變化。 這樣的研究已經(jīng)有很多， 包括不同 時(shí)間點(diǎn)或者不同空間取樣來分析腫瘤異質(zhì)性的，比如： Marusyk A, Almendro V, Polyak K. Intra-tumour heterogeneity: a looking glass for cancer?J. Nature Reviews Cancer, 2012, 12(5): 323-334.Zhang J, Fujimoto J, Zhang J, et al. Intratumor heteroge

3、neity in localized lung adenocarcinomas delineated by multiregion sequencingJ. Science, 2014, 346(6206): 256-259.Shi Y J, Tsang J, Ni YB, et al. Intratumoral Heterogeneity in Breast Cancer: A Comparison of Primary and Metastatic Breast CancersJ. Oncologist, 2017, 22(4).Hao J J, Lin D C, Dinh H Q, et

4、 al. Spatial intratumoral heterogeneity and temporal clonal evolution in esophageal squamous cell carcinomaJ.Nature genetics, 2016, 48(12): 1500-1507.下面就拿 Naturegenetics, 2016 的關(guān)于 ESCC 不同空間取樣探究腫瘤異質(zhì)性 來詳細(xì)解讀。文章解讀食管癌， Esophageal squamous cell carcinoma (ESCC) ，作者選擇了 13 個(gè) ESCC 病人的 51 處 腫瘤組織進(jìn)行 WES 測(cè)序分析，同時(shí)

5、也測(cè)了他們的正常組織 做對(duì)照。平均測(cè)序深度 150X ，測(cè)序策略都是 150PE ?？偣?找到了涉及了 1427 個(gè)基因的 1610 個(gè)非沉默突變， 還有 568 個(gè)沉默突變。本文主要關(guān)注點(diǎn)就是： spatial intratumoral heterogeneity (ITH) and temporal clonal evolutionarysilent mutations (沉默突變)：即同義突變， 突變雖然替換了堿基， 但氨基酸順序未變，保持野生型的功能。有趣的是其中兩個(gè) 病人用的是 BGI 的 CG 測(cè)序平臺(tái)，而其余的都是用 Agilent SureSelect Human All Ex

6、on v4 (51 Mb) kit捕獲外顯子序列，用 Illumina HiSeq 4000 進(jìn)行 PE150 的測(cè)序，當(dāng)然，還 是在 BGI 公司測(cè)的。他們選擇的分析流程是最經(jīng)典的 BWA-MEM 比對(duì)到 hg19 ，然后走 GATK best practices 得 到每個(gè)測(cè)序樣本的 bam 文件。至于找 somatic mutation 步驟，選取的是 varscan ，注釋用的是 ANNOVAR 。對(duì)于得到腫瘤特異性的變異之后做的高級(jí)分析包括： Phylogenetic tree constructionCancer cell fraction analysisIdentificati

7、on of putative driver mutationsMutational signature analysisDNA methylation analysis 甲基化數(shù)據(jù) (Illumina HumanMethylation450 BeadChip) 都上傳到了 GEO 里面， GSE79366, 對(duì)應(yīng)的 NGS 序列也在 SRA 數(shù)據(jù)庫中 SRP072112 ，我們可以獲取到數(shù)據(jù)。我隨便選了一個(gè)人的 四個(gè)全外顯子數(shù)據(jù)，測(cè)序策略都是 PE150 ， SRA 數(shù)據(jù)庫的 ID 分別是 :ESCC13-T4 SRR3270888ESCC13-T3 SRR3270887ESCC13-T2 S

8、RR3270886ESCC13-T1 SRR3270885ESCC13-N SRR3270884 從 SRA 數(shù)據(jù)庫 下載并轉(zhuǎn)換為 fastq 測(cè)序數(shù)據(jù)文件把上面的描述文本存為文 件 下載腳本如下： | cut -f 2|while read iddo echo $idwget -c dy/sra/SRP/SRP072/SRP072112/$id/$id.sradone 轉(zhuǎn)換腳 本如下 | -gzip -split-3 -A $array0 $array1.sradone 得到的 sra 和 fastq 文件如下： 5.7G Sep 19 17:2419 17:1819 16:1619 15:

9、5919 15:5020 10:3820 10:38 20 11:32 20 11:32 20 06:44 20 06:44 20 03:30 20 03:30 20 00:42 20 00:42大家看上面文件的日期，可以很明 顯的看到我下載 sra 文件，以及把 sra 文件用 fastq-dump 轉(zhuǎn) 為 fastq 文件所消耗的時(shí)間， 基本上下載只需要 15 分鐘， 但 是轉(zhuǎn)換居然耗時(shí) 3 小時(shí)，所以并行很有必要哈。所以一晚上 才轉(zhuǎn)換了 T2,T3,T4 這 3 個(gè)樣本，我嫌棄它太慢了，就把 T1 和N兩個(gè)樣本并行了。簡(jiǎn)單的走一下fastqc+multiqc 看看數(shù) 據(jù)質(zhì)量，一般都會(huì)很不

10、錯(cuò)的，這個(gè)數(shù)據(jù)也不例外。 ls *.gz |xargs /biosoft/fastqc/FastQC/fastqc -o ./ -t 5 然后走 WES的標(biāo)準(zhǔn) SNP-calling 流程選用的是經(jīng)典的 GATK best practice 的流程，整個(gè)項(xiàng)目最后耗費(fèi)空間約 500G ，代碼如 下 :#!/bin/bash#SBATCH -job-name wes_tumor_human#SBATCH -partition FHS_NORMAL#SBATCH -nodes 1#SBATCH -tasks-per-node5#SBATCH -mem40G #SBATCH -outputwes.%j

11、.out#SBATCH-errorwes.%j.err#SBATCH -mail-typeFAIL#SBATCH -(using htslib # # Step 1 : Alignment # #start=$(date +%s.%N)echo bwa datebwa mem -t 5 -M -R RGtID:$sampletSM:$sampletLB:WEStPL:Illumina $INDEX $fq1 $fq2 1$sample.sam 2/dev/null echo bwa datedur=$(echo $(date +%s.%N) - $start | bc)printf Execu

12、tion time for BWA : %.6f seconds $durecho # # Step 2: Sort and Index # #start=$(date +%s.%N)echo-Xmx40g -jar $PICARD SortSam SORT_ORDER=coordinate INPUT=$sample.sam OUTPUT=$sample.bamsamtools index $sample.bamecho SortSam datedur=$(echo $(date +%s.%N) - $start | bc)printf Execution time for SortSam

13、: %.6f seconds $durecho rm $sample.sam# # Step 3: Basic Statistics # #start=$(date +%s.%N)echo stats datesamtools flagstat $sample.bam $sample.alignment.flagstatsamtools stats $sample.bam $sample.alignment.statecho plot-bamstats -p $sample_QC $sample.alignment.statecho stats datedur=$(echo $(date +%

14、s.%N) - $start | bc)printf Execution time for Basic Statistics : %.6f seconds $durecho # # Step 4: multiple filtering for bam files # #MarkDuplicates#start=$(date+%s.%N)echo MarkDuplicates datejavaMarkDuplicates INPUT=$sample.bamOUTPUT=$sample_marked.bamMETRICS_FILE=$sample.metricsecho MarkDuplicate

15、s datedur=$(echo $(date +%s.%N) - $start | bc)printf Execution time for MarkDuplicates : %.6f seconds $durecho rm $sample.bam-Xmx40g -jar $PICARD FixMateInformation INPUT=$sample_marked.bam OUTPUT=$sample_marked_fixed.bam SO=coordinatesamtools index $sample_marked_fixed.bamecho FixMateInfo datedur=$

16、(echo $(date +%s.%N) - $start | bc)printf Execution time for FixMateInfo : %.6f seconds $durecho rm $sample_marked.bam# # Step 5: gatk process bam files# # SplitNCigar #start=$(date+%s.%N)echo SplitNCigar datejavaReassignOneMappingQuality -RMQF 255 -RMQT 60 -U ALLOW_N_CIGAR_READS#-fix_misencoded_qua

17、lity_scores#-fix_misencoded_quality_scores only if phred 64 echo SplitNCigar datedur=$(echo $(date +%s.%N) - $GATK -T RealignerTargetCreator -I $sample_marked_fixed_split.bam -R $GENOME -o $sample_ervals -known $Mills_indels -known $KG_indels -nt 5echo RealignerTargetCreator datedur=$(echo

18、 $(date +%s.%N) - $start | bc)printf Execution time for RealignerTargetCreator : %.6f seconds$durecho #IndelRealigner#start=$(date$sample_recal.bamecho PrintReads datedur=$(echo $(date +%s.%N) - $start | bc)printf Execution time for PrintReads : %.6f seconds $durecho rm $sample_realigned.bamchmod ug

19、a=r $sample_recal.bam# # Step 6: gatk call snp/indel# # start=$(date- -o $sample_raw_indels.vcf# :# -filterExpression QD60.0 | MQ-filterName my_snp_filter -V$sample_raw_snps.vcf-o $sample_filtered_snps.vcf-T SelectVariants -R $GENOME-excludeFiltered -VVariantFiltration -R $GENOME-filterExpression QD

20、200.0 | ReadPosRankSum-filterName my_indel_filter -V $sample_raw_indels.vcf -ojava-Xmx40g -jar $GATK -T -excludeFiltered -V -o把上面的代碼保存 為 humantumorwes.sh 文件，運(yùn)行即可： 比對(duì)成功后得到的 sam/bam 文件如下； 38G Sep 20 15:15 ESCC13-N.sam52G Sep 20 16:35 ESCC13-T1.sam43G Sep 20 11:34 ESCC13-T2.sam38G Sep 20 10:39 ESCC13-T

21、3.sam40G Sep 20 12:34 ESCC13-T4.sam# 本來只需要一天就可以完成的，因?yàn)橹虚g不小心誤刪了，所 以又花費(fèi)了一天重新運(yùn)行。 17G Sep 22 00:01 ESCC13-N_recal.bam22G Sep 22 02:09ESCC13-T1_recal.bam21G Sep 22 01:59 ESCC13-T2_recal.bam16G Sep 21 23:29 ESCC13-T3_recal.bam17G Sep 22 00:32 ESCC13-T4_recal.bam 用 GATK 做 Snp-calling 結(jié)束后得到 的 germline 的 vcf

22、如下： 18M Sep 22 05:14 ESCC13-N_filtered_indels.vcf 15M Sep 22 05:14 22 05:13 22 05:13 ESCC13-N_filtered_snps.vcf 17M Sep 22 05:11 ESCC13-N_raw_indels.vcf163M Sep 22 05:10 22 05:11 ESCC13-N_raw_snps.vcf261M Sep 22 08:45 22 08:27 ESCC13-T2_filtered_indels.vcf 21M Sep 22 08:27 22 08:26 22 08:25 ESCC13-T

23、2_filtered_snps.vcf 24M Sep 22 08:23 ESCC13-T2_raw_indels.vcf240M Sep 22 08:22 22 08:23 ESCC13-T2_raw_snps.vcf 19M Sep 22 05:03 ESCC13-T3_filtered_indels.vcf 16M Sep 22 05:04 22 05:0222 05:02 ESCC13-T3_filtered_snps.vcf 18M Sep 22 05:01 ESCC13-T3_raw_indels.vcf169M Sep 22 04:59 22 05:00 ESCC13-T3_raw_snps.vcf 17M Sep 22 06:02 ESCC13-T4_filtered_indels.vcf 14M Sep 22 06:03 22 06:01 22 06:01 ESCC13-T4_filtered_snps.vcf 17M Sep 22 06:00 ESCC13-T4_raw_indels.vcf156M Sep 22 05:59 22 05:59 ESCC13-T4_raw_snps.vcf 可以看到同一個(gè)病人的不同部位 的測(cè)試數(shù)據(jù)從 VCF 文件大