解码生命 守护健康

实战:同一肿瘤内 基因突变异质性简单探究

2017-09-19 08:03:19生信之光

大家好!我是阿尔的太阳,这次是实战篇,我们将在SRA数据库上下载公共数据,利用服务器对其进行分析,并对感兴趣的问题作初步的探究


今天的主题是同一肿瘤内基因突变异质性简单探究

探究在同一个时间点的同一个肿瘤病人的同一个肿瘤组织的

不同部位的WES测序,看看均一性如何


均一性越低说明,肿瘤内部的异质性越高


实战三步走:

数据获取 进行分析 初步结论 


1.数据获取

数据来源是 Nat Genet. 2016 December  的一篇文章

主题是鼻咽癌的不同病人的肿瘤内的异质性研究


在SRA数据库中我们可以获取到数据 


我随便选了一个人的四个全外显子数据

分别是


SRR3270884-正常组织

SRR3270885-肿瘤组织1

SRR3270886-肿瘤组织2

SRR3270887-肿瘤组织3

SRR3270888-肿瘤组织4


感兴趣可以到SRA数据库搜索SRR ID 进行查看 


可以看到其平台为 Illumina Hiseq 4000 [ X TEN ] 


模式为PE150


下载的代码是download_SRA.sh 

SRR_id=$1

target_dir=$2


#ascp=/home/ytan/.aspera/connect/bin/ascp

#ascpid=/home/ytan/.aspera/connect/etc/asperaweb_id_dsa.openssh

ascp=/home/glad_zarl/.aspera/connect/bin/ascp

ascpid=/home/glad_zarl/.aspera/connect/etc/asperaweb_id_dsa.openssh

NCBI_SRAlink=anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/SRR


$ascp -i $ascpid \

-k 1 -T -l200m ${NCBI_SRAlink}/${SRR_id:0:6}/${SRR_id}/${SRR_id}.sra \

$target_dir > ${target_dir}/${SRR_id}.log 


解压的代码是:sra_uncom.sh

#the src's mission is sra_uncom BGI seq

dir=$1

tardir=$2


#sratoolkit

sratoolkit=/thinker/net/biosoft/sratoolkit/sratoolkit.2.8.1-1-centos_linux64/bin/fastq-dump


echo "$dir"

cd $dir


#do it

for i in *.sra

do

echo "begin to umcom $i"

time \

$sratoolkit \

--split-3 \

$i \

--outdir $tardir > ${tardir}/${i}uncom.log


echo "finish to umcom $i"

done


就得到了 raw fastq


初始数据大小

5.7G 9月  16 23:03 SRR3270884.sra

7.7G 9月  17 12:58 SRR3270885.sra

6.8G 9月  17 12:56 SRR3270886.sra

5.7G 9月  16 23:47 SRR3270887.sra

6.0G 9月  17 14:02 SRR3270888.sra


  2.进行分析 数据分析:

2.1

首先是fastq过滤掉低质量的reads,在这里我们会用到海普洛斯开发的afterqc软件进行自动化过滤


过滤的代码是:After_qc.sh 

fq=$1

o=$2


python /thinker/net/zarl-2/STUDIES/WES/sra_fq/AfterQC-master/after.py \

--no_correction \

-1 ${fq}_1.fastq \

-2 ${fq}_2.fastq \

-g $o \

-b $o \

-r $o


fastqQC-itools.sh

#the src's mission is to fastqc

dir=$1

tardir=$2


#itools

itools=/thinker/net/biosoft/REseqtools/Reseqtools/iTools_Code/iTools

echo "$dir"

cd $dir


#do it

#for i in *.fastq.gz

for i in *.fastq

do

echo "begin to itools $i"

time $itools Fqtools stat -InFq $i -OutStat ${tardir}/${i}_itools.txt

echo "finish to itools $i"


done


过滤之后的数据大小

17G 9月  17 23:16 N_1.fastq

17G 9月  17 23:16 N_2.fastq

23G 9月  18 00:36 T1_1.fastq

23G 9月  18 00:36 T1_2.fastq

19G 9月  17 23:53 T2_1.fastq

19G 9月  17 23:53 T2_2.fastq

17G 9月  17 23:30 T3_1.fastq

17G 9月  17 23:30 T3_2.fastq

17G 9月  17 23:27 T4_1.fastq

17G 9月  17 23:27 T4_2.fastq


过滤之后的数据质量

可以看到数据质量非常好

可以用于后续的分析步骤


2.2

比对到参考基因组 转换为bam 并sort 好 并去除 pcr  重复

基于经验来说  基于探针捕获的测序,PCR重复会有相当一部分

但是WGS的PCR重复会很少 尤其是 PCR-FREE 的测序

比对和排序:imported-fq-to-raw-bam.sh

name=$1

outdir=$2

#get the date of the task begining

start_time=$(date +%H:%M:%S)

echo "Begin at $start_time"


#use the bwa mem to mapping reads to the ref genome 

#ref genome format is hg19

echo "1.mapping the NGS reads to the ref genome"

glad bwa mem -t 4 -M -R "@RG\tID:081-300m\tSM:${name}\tLB:PE\tPL:Illumina" \

  /thinker/dstore/world/gene/sapiens/reference/hututa/glad-ref-r1/v2/human_g1k_v37 \

  ${name}_1.fastq \

  ${name}_2.fastq \

  > $outdir/${name}.sam


echo "2.process the bam file : sort & refine & index"

#sort the bam and index using the samtools

glad samtools sort -@ 8 $outdir/${name}.sam > $outdir/${name}-sorted.bam

glad samtools index $outdir/${name}-sorted.bam


去重复:rmdup_bam.sh

#the src's mission is filter bam file rmdup

# exit 1 if any err

set -o errexit


#mission is rmdup

name=$1

dir=$2

tar=$3

nohup java -jar /thinker/net/biosoft/picardtools/picard-tools-1.119/MarkDuplicates.jar \

I=${dir}${name}.bam O=${tar}${name}rmdup.bam \

M=${tar}${name}marked_dup_metrics.txt \

REMOVE_DUPLICATES=true \

ASSUME_SORTED=true > ${tar}${name}mark.log &

echo "$i rmdup end"


2.3

使用speedseq  somatic call vars 

https://github.com/hall-lab/speedseq 

安装和使用十分简单

/thinker/net/zarl-2/zarlsoft/speedseq/bin/speedseq somatic \

-t 24 /thinker/dstore/world/gene/sapiens/reference/hututa/glad-ref-r1/v2/human_g1k_v37.fa \

N-sortedrmdup.bam T4-sortedrmdup.bam > T4.log


2.4

使用annovar 对somatic mutation 进行基本的注释

#cover format

${annovar}/convert2annovar.pl \

 -format vcf4old $vcf > ${dir}/${name}_anno/${name}.annovar


${annovar}/annotate_variation.pl -buildver hg19 --geneanno --outfile ${dir}/${name}_anno/${name}.anno ${dir}/${name}_anno/${name}.annovar $db

echo "ok!"


2.5

使用shell 抓取基因名称

cat T1-sortedrmdup.bam.vcf.anno.exonic_variant_function | cut -f 3 | tr ":" "\t" | cut -f 1 | sort | uniq > T1.list


2.6

懒得写代码了

使用在线VEEN图检测均一性 


http://bioinformatics.psb.ugent.be/cgi-bin/liste/Venn/calculate_venn.htpl 


3.初步结论


可以看到T1体细胞突变涉及的基因是2546个,然后以此类推

可以看到,这四块组织涉及的突变的基因数量是不同的

T4最多,T2最少


List namesnumber of elementsnumber of unique elements
T1.list25462546
T2.list25352535
T3.list24052405
T4.list33133313
Overall number of unique elements5280

可以看到 T1-T3 涉及到的突变基因一致性比较高 但是第四块组织却涉及到了更加多的突变的基因数量 


即使同一个人的同一个时间点的同一块组织的几个不同的位置

基因突变也是不一样的 

可见肿瘤是一个高度异质性的,组织水平的疾病

也可以看到单时间点  单块组织测序的局限性 

其只能反映一个局部区域的基因突变特征