实战:同一肿瘤内 基因突变异质性简单探究
大家好!我是阿尔的太阳,这次是实战篇,我们将在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 names | number of elements | number of unique elements |
T1.list | 2546 | 2546 |
T2.list | 2535 | 2535 |
T3.list | 2405 | 2405 |
T4.list | 3313 | 3313 |
Overall number of unique elements | 5280 |
可以看到 T1-T3 涉及到的突变基因一致性比较高 但是第四块组织却涉及到了更加多的突变的基因数量
即使同一个人的同一个时间点的同一块组织的几个不同的位置
基因突变也是不一样的
可见肿瘤是一个高度异质性的,组织水平的疾病
也可以看到单时间点 单块组织测序的局限性
其只能反映一个局部区域的基因突变特征