肿瘤外显子数据分析指南 复现
文章目录
- 一、读文献并下载数据
- 1.1、试图使用sratoolkit下载，虽然总是失败，但原理上应该没错
- 1.2、获取SRP070662的下载链接
- 二、质控与去接头
- 2.1 fastqc
- 2.2 质控和去接头
- 3、比对与bam文件格式
- 三、比对结果的质控
- 四、GATK流程和找变异
- 4.1标记或去除重复序列
- 4.2 碱基质量重校正BQSR
- 三、单样本找Germline SNPS+INDELs
- 四、mutect找Somatic mutations
一、读文献并下载数据
1.1、试图使用sratoolkit下载，虽然总是失败，但原理上应该没错
#安装sratoolkit,从NCBI主页获取软件包下载地址
1.1.1、下载sratoolkit软件包
wget -P ~/SCS2021/wes_cancer/biosoft https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.11.0/sratoolkit.2.11.0-centos_linux64.tar.gz
1.1.2、解压压缩包
#用tar命令的-C参数，设置解压文件保存路径在 ~/Scs2021/wes_cancer/biosoft
tar -zvxf sratoolkit.2.11.0-centos_linux64.tar.gz -C .
1.1.3、测试安装是否成功
~/SCS2021/wes_cancer/biosoft/sratoolkit.2.11.0-centos_linux64/bin/fastq-dump -h
#4、将sratoolkit安装文件路径加入环境变量
echo ‘export PATH=~/SCS2021/wes_cancer/biosoft/sratoolkit.2.11.0-centos_linux64/bin:$PATH’ >> ~/.bashrc
source ~/.bashrc
1.2、获取SRP070662的下载链接
#1.2.1 wget方法下载（本人采用，因为听说wget更加稳定，缺点是很慢）
去ENA主要搜索SRP070662，下载case1
将SRR号都写在~/SCS2021/wes_cancer/project/SRR_Acc_List.txt中，cat后展示如下
#批量获得SRR_1代码
cat SRR_Acc_List.txt|while read id ;do t1=${id:0:6};t2=${id:9:10};t3=${id:0:10};echo "ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/fastq/$t1/00$t2/$t3/${t3}_1.fastq.gz ." ;done
#批量获得SRR_2代码
cat SRR_Acc_List.txt|while read id ;do t1=${id:0:6};t2=${id:9:10};t3=${id:0:10};echo "ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/fastq/$t1/00$t2/$t3/${t3}_2.fastq.gz ." ;done
#写入fq_2.sh中
给下载好的fastq文件重命名
最终结果
#1.2.2 aspera下载（尝试，但总是断，于是放弃）
#去ENA主要搜索SRP070662，选择aspera-sra数据，下面会出来相应下载网址，其余步骤参考1.2.1，不具体写了，以下代码供参考（不完整）
#得出批量下载代码
`cat SRR_Acc_List.txt|while read id
do
t1=${id:0:6};t2=${id:9:10};t3=${id:0:10}
echo "ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/fastq/$t1/00$t2/$t3/$t3.fastq.gz ."
done#得出批量转移代码
cat SRR_Acc_List.txt|while read id
do t1=${id:0:6};t2=${id:9:10};t3=${id:0:10}
echo "mv ${t3}_1.fastq.gz ${t3}_exome_sequencing_of_case1_germline_1.fastq.gz"
done#写入一个脚本中运行
cd ~/SCS2021/wes_cancer/project/1.raw
touch fq_download.sh
vi fq_download.sh #放入上述批量下载代码和转移代码
#!/usr/bin/env bash
nohup bash fq_download.sh &3.迅雷下载
注：当用wget下载不能确定是否下载完整的时候，可以使用迅雷看差不多有多少个G
**
二、质控与去接头
**
2.1 fastqc
在project文件夹下创建config，进行批量操作获得相应代码
cat config | while read id
do
fastqc --outdir ./3.qc/raw_qc --threads 16 ./1.raw_fq/${id}*.fastq.gz >> ./3.qc/raw_qc/${id}_fastqc.log 2>&1
donemultiqc ./3.qc/raw_qc/*.zip -o ./3.qc/raw_qc/multiqc
#经过fastqc处理后，每个fastq文件会生成一个html文件和一个zip压缩文件，其中html文件是指控报告，可以下载到本地用浏览器打开查看，而zip文件里面是html文件图表对应的数据。因为样本数较多，我们就用multiqc软件对所有的样本的结果进行整合，方便一起查看。
结果：
2.2 质控和去接头
touch trim_galore.sh
vi trim_galore.sh
cat SRR_Acc_List.txt | while read id
do
fq1=./1.raw_fq/${id}_1.fastq.gz
fq2=./1.raw_fq/$[id}_2.fastq.gz
trim_galore --paired -q 28 --phred33 --length 30 --stringency 3 --gzip --cores 8 -o ./2.clean_fq $fq1 $fq2 >> ./2.clean_fq/${id}_trim.log 2>&1
donenohup bash trim_galore.sh &
3、比对与bam文件格式
3.1 构建索引
conda activate wes
cd ~/wes_cancer/data
gunzip Homo_sapiens_assembly38.fasta.gz
time bwa index -a bwtsw -p gatk_hg38 ~/wes_cancer/data/Homo_sapiens_assembly38.fasta
cd ~/wes_cancer/project
3.2 比对
## bwa.sh
cat SRR_Acc_List.txt | while read id
do
echo "start bwa for ${id}"
fq1=./2.clean_fq/${id}_1_val_1.fq.gz
fq2=./2.clean_fq/$[id}_2_val_2.fq.gz
bwa mem -M -t 24 -R "@RG\tID:${id}\tSM:${id}\tLB:WXS\tPL:Illumina" ~/SCS2021/wes_cancer/data/gatk_hg38 ${fq1} ${fq2} | samtools sort -@ 24 -m 1G -o ./4.align/${id}.bam
echo "end bwa for ${id}"
done
#bwa比对后先产生sam文件，但是因为sam文件太大，占用空间，要转换成二级制的bam文件，两者记录的信息是一样的。因此，通常都是把比对的结果通过管道符|传给samtools进行排序后直接生成bam文件，排序的依据是bam文件中reads比对到参考基因组的染色体和坐标。生成的bam文件大小为3~13G不等。
最终结果
3.3 提取小bam
touch smallbam.sh
vi smallbam.sh
cd 4.align
#!/usr/bin
cat SRR_Acc_List.txt|while read id
do
samtools index ${id}.bam
samtools view -h ${id}.bam chr17|samtools view -Sb -> ${id}_small.bam
samtools index ${id}_small.bam
done载入IGV
三、比对结果的质控
##得到了bam文件之后，也要进行质控，可以对测试数据有更深入的了解。
cd SCS2021/wes_cancer/project
touch stats.sh
cat SRR_Acc_List.txt | while read id
dobam=./4.align/${id}.bamsamtools stats -@ 16 --reference ~/SCS2021/wes_cancer/data/Homo_sapiens_assembly38.fasta ${bam} > ./4.align/stats/${id}.stat##plot-bamstats一直报错，不知为何plot-bamstats -p ./4.align/stats/${id} ./4.align/stats/${id}.stat
done
这一步质控后也会生成html网页报告，下载到本地查看，可以看到碱基分布信息和统计信息。
##使用qualimap软件来查看基因组覆盖深度等信息
touch qualimap.sh
vi qualimap.sh
cat SRR_Acc_List.txt | while read id
do
qualimap bamqc --java-mem-size=10G -gff ~/SCS2021/wes_cancer/data/hg38.exon.bed -nr 100000 -nw 500 -nt 16 -bam ./4.align/${id}.bam -outdir ./4.align/qualimap/${id}
done
四、GATK流程和找变异
我们已经从原始的fastq文件经trim过滤后，再用bwa比对到参考基因组得到了bam文件，所以接下来走Mark Duplicates和BQSR
4.1标记或去除重复序列
在拿到的测序结果中，有一部分序列是重复的，这是由于上机前xuyaojinxingPCR扩增，如果序列扩增次数不同，就会导致PCR重复#标记或去除重复序列
vi gatk.sh
GATK=~SCS2021/wes_cancer/biosoft/gatk-4.2.0.0/gatkcat SRR_Acc_List.txt | while read id
doBAM=./4.align/${id}.bamif [ ! -f ./5.gatk/ok.${id}_marked.status ]thenecho "start MarkDuplicates for ${id}" `date`$GATK --java-options "-Xmx20G -Djava.io.tmpdir=./" MarkDuplicates \-I ${BAM} \--REMOVE_DUPLICATES=true \-O ./5.gatk/${id}_marked.bam \-M ./5.gatk/${id}.metrics \1>./5.gatk/${id}_log.mark 2>&1 if [ $? -eq 0 ]thentouch ./5.gatk/ok.${id}_marked.statusfiecho "end MarkDuplicates for ${id}" `date`samtools index -@ 16 -m 4G -b ./5.gatk/${id}_marked.bam ./5.gatk/${id}_marked.baifi
donevi samtools.sh##if判断语句判断是否已经存在生成的./5.gatk/ok.${id}_marked.status文件，如果存在，就不再执行，如果不存在就继续往下走。如果在shell命令执行后，会有一个返回值$?，如果返回值为0，只上一步命令成功执行，否则执行失败。在这里，如果命令成功执行，就创建一个空文件 ok.${id}_marked.status
4.2 碱基质量重校正BQSR
###BaseRecalibrator 根据原始bam文件中碱基质量值计算出系统误差的分布.这里计算了所有需要进行重校正的reads和特征值，然后把这些信息输出为一份校准表文件（recal.table）。--know-sites参数需要输入已知且可靠的变异位点
###ApplyBQSR 根之前计算的模型对碱基质量进行校正。这一步利用得到的校准表文件(recal.table)重新调整原来BAM文件中的碱基质量值，并使用这个新的质量值重新输出一份新的BAM文件
vi bqsr.sh
GATK=~/SCS2021/wes_cancer/biosoft/gatk-4.2.0.0/gatk
snp=~/SCS2021/wes_cancer/data/dbsnp_146.hg38.vcf.gz
indel=~/SCS2021/wes_cancer/data/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
ref=~/SCS2021/wes_cancer/data/Homo_sapiens_assembly38.fasta
cat SRR_Acc_List.txt | while read id
doif [ ! -f ./5.gatk/${id}_bqsr.bam ]thenecho "start BQSR for ${id}" `date`$GATK --java-options "-Xmx20G -Djava.io.tmpdir=./" BaseRecalibrator \-R $ref \-I ./5.gatk/${id}_marked.bam \--known-sites ${snp} \--known-sites ${indel} \-O ./5.gatk/${id}_recal.table \1>./5.gatk/${id}_log.recal 2>&1 fi
donevi bqsr_2.sh
GATK=~/SCS2021/wes_cancer/biosoft/gatk-4.2.0.0/gatk
snp=~/SCS2021/wes_cancer/data/dbsnp_146.hg38.vcf.gz
indel=~/SCS2021/wes_cancer/data/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
ref=~/SCS2021/wes_cancer/data/Homo_sapiens_assembly38.fasta
cat SRR_Acc_List.txt | while read id
do$GATK --java-options "-Xmx20G -Djava.io.tmpdir=./" ApplyBQSR \-R $ref \-I ./5.gatk/${id}_marked.bam \-bqsr ./5.gatk/${id}_recal.table \-O ./5.gatk/${id}_bqsr.bam \1>./5.gatk/${id}_log.ApplyBQSR 2>&1 echo "end BQSR for ${id}" `date`
done
三、单样本找Germline SNPS+INDELs
GATK=~/SCS2021/wes_cancer/biosoft/gatk-4.2.0.0/gatk
snp=~/SCS2021/wes_cancer/data/dbsnp_146.hg38.vcf.gz
indel=~/SCS2021/wes_cancer/data/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
ref=~/SCS2021/wes_cancer/data/Homo_sapiens_assembly38.fasta
bed=~/SCS2021/wes_cancer/data/hg38.exon.bedcat SRR_Acc_List.txt | while read id
doecho "start HC for ${id}" `date`$GATK --java-options "-Xmx20G -Djava.io.tmpdir=./" HaplotypeCaller -ERC GVCF \-R ${ref} \-I ./5.gatk/${id}_bqsr.bam \--dbsnp ${snp} \-L ${bed} \-O ./5.gatk/${id}_raw.vcf \1>./5.gatk/${id}_log.HC 2>&1echo "end HC for ${id}" `date`donecd ./5.gatk/gvcf
for chr in chr{1..22} chrX chrY chrM
dotime $GATK --java-options "-Xmx20G -Djava.io.tmpdir=./" GenomicsDBImport \
-R ${ref} \
$(ls ./*raw.vcf | awk '{print "-V "$0" "}') \
-L ${chr} \
--genomicsdb-workspace-path gvcfs_${chr}.dbtime $GATK --java-options "-Xmx20G -Djava.io.tmpdir=./" GenotypeGVCFs \
-R ${ref} \
-V gendb://gvcfs_${chr}.db \
-O gvcfs_${chr}.vcfdone$GATK --java-options "-Xmx20G -Djava.io.tmpdir=./" GatherVcfs \
$(for i in {1..22} X Y M;do echo "-I gvcfs_chr${i}.vcf" ;done) \
-O merge.vcf
四、mutect找Somatic mutations
###前面单样本走GATK HaplotypeCaller流程得到vcf文件主要包含了germline mutations 信息。但是我们这里分析的是肿瘤样本，包含成对的tumor和normal样本，更加关心的是tumor样本中所包含的somatic muattions信息。所以我们应该走GATK的somatic mutations流程###我们这里有成对的normal和tumor样本，需要把配置文件config2做成两列
$ cat config2
case1_germline case1_biorep_A_techrep_1
case1_germline case1_biorep_B
case1_germline case1_biorep_C
case1_germline case1_techrep_2
未完待续…
肿瘤外显子数据分析指南 复现相关推荐
- TMF大数据分析指南 Unleashing Business Value in Big Data（一）
大数据分析指南 TMF Frameworx最佳实践 Unleashing Business Value in Big Data 前言 此文节选自TMF Big Data Analytics Guide ...
- 【数据分析】Python数据分析指南(全)
前言 数据分析是通过明确分析目的,梳理并确定分析逻辑,针对性的收集.整理数据,并采用统计.挖掘技术分析,提取有用信息和展示结论的过程,是数据科学领域的核心技能. 本文从数据分析常用逻辑框架及技术方法出 ...
- 恒瑞医药卡瑞利珠单抗七大瘤种获2022版中国临床肿瘤学会诊疗指南19项推荐
上海2022年4月30日 /美通社/ -- 近日,<中国临床肿瘤学会常见恶性肿瘤诊疗指南(2022版)>(简称"CSCO指南")完成更新并已发布.基于循证医学证据.汇聚 ...
- 超级详细的 Python 数据分析指南
来源 | 算法进阶 头图 | 下载于视觉中国 前言: 数据分析是通过明确分析目的,梳理并确定分析逻辑,针对性的收集.整理数据,并采用统计.挖掘技术分析,提取有用信息和展示结论的过程,是数据科学领域的核 ...
- SEER 肿瘤数据库使用指南
SEER 数据库是由美国国立癌症研究所于 1973 年建立,是美国常用的癌症数据库,里面包括各式各样的肿瘤类型,如肺癌.乳腺癌.胃癌.结直肠癌.前列腺癌等等. 主要提供了各式各样的临床资料,如性别.年 ...
- Python数据分析指南
公众号:尤而小屋 作者:Peter 编辑:Peter 大家,我是Peter~ 关注尤而小屋的朋友都知道,Peter一直都在写关于数据分析的文章. 有读者朋友后台反映:Peter能不能将公众号的文章进行 ...
- 实用流量数据分析指南
本文转自公众号『数据管道』,详情请扫码关注该公众号: 前戏 粽子节了嘛,突然想吃粽子了,咋办,买粽子呗!现在情景转换一下,假设你是某饮食网的数据分析师,现在某粽子界大亨想拿钱砸你老板(打广告). 老板 ...
- 肿瘤筛查指南--依据薄世宁的得到讲座总结
肿瘤筛查123: 1种细菌:幽门螺旋杆菌 2种病毒: HPV 病毒(女性).乙肝病毒: 3种必要检查:男性要查肺部低剂量 CT .胃肠镜.前列腺超声+ PSA :女性要查肺部低剂量 CT .胃肠镜.乳 ...
- anaconda pycharm_搭建 Python 高效开发环境： Pycharm + Anaconda
介绍 先来介绍下两位主角: Pycharm:目前一款主流的 Python 集成开发环境,它带有一整套帮助我们在Python开发时提高效率的工具,比如调试.语法高亮.Project管理.代码跳转.智能提 ...
- python write 写多行_如何用 Python 执行单行命令
一般来说,面对日常处理的一些小任务,直接用 sed,grep 之类的就可以搞定,更复杂一点的就会考虑 awk 或者用一些现成的轮子,要是 awk 搞不定我就只好用 Python 了.但有些时候,我仅仅 ...
最新文章
- Android开发之动态库调用
- 驱动学习之LED驱动框架
- Android数据适配器Adapter简介
- 探讨磷酸铁锂电池在UPS的应用
- SQL查询语句精华文章(转)
- 历史数据如何处理_数据库表数据量大读写缓慢如何优化（1）【冷热分离】
- 一行c语言代码,打钩的一行c语言代码解释一下,谢谢,详细解释绝对最佳
- 收藏 | Python必备技能之 25个Matplotlib常用代码！
- Elasticsearch 读时分词、写时分词
- 51单片机——LCD12864
- Excel如何将多个sheet导出到PDF？
- 顶顶通软电话介绍-一个网络电话客户端（SIP软电话）
- H5 在iPhone真机上调试H5页面
- GMap.NET 使用教程【1】
- go语言黑帽子学习3
- vue中使用ECharts实现折线图和饼图
- 任正非谈“咖啡杯”文化
- Windows下通过注册表修改某个类型文件的默认打开方式和文件图标
- #include int inc(int a) { return(++a); } int multi(int*a,int*b,int*c) { return(*c=*a**b); }
- 【机器学习】课程笔记16_大规模机器学习(Large Scale Machine Learning)