2019-12-21

12.21Wes -seq笔记总结1

上午 wes-seq

之前跑了rna-seq的上游总结了流程，现在想把上周的东西整合归纳一下，在跑一遍

###建各级目录
mkdir data biosoft project 
cd project
mkdir -p 1.raw_fq  2.qc/{raw_qc,clean_qc} 3.clean_fq  4.align/{qualimap,featurecount,stat,flagstat}  5.gatk/gvcf   6.mutect2   7.annovar  8.cnvkit   9.facet #建立各级目录
tree -L 1 #查看目录

软件准备

安装conda

介绍过很多次了,这里不说啦

###安装conda
wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh 
bash Miniconda3-latest-Linux-x86_64.sh 
source ~/.bashrc#永久激活
# 下面这三行配置官网的channel地址
conda config --add channels r 
conda config --add channels conda-forge 
conda config --add channels bioconda
# 下面这四行配置清华大学的bioconda的channel地址
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/
conda config --set show_channel_urls yes

创建wes小环境

conda  create -n wes python=2
conda activate wes#激活小环境

## biosoft
conda install -y fastqc trim-galore samtools bwa vcftools bedtools cnvkit subread qualimap #批量下载软件
##记得调用help，可以调用就ok

## GATK
wget -c https://github.com/broadinstitute/gatk/releases/download/4.1.4.0/gatk-4.1.4.0.zip
#直接下载就可以使用，不需要安装

用到的数据

GATK只能用GATK官网的参考基因组，所以只能从其官网下载

## GATK所用到的数据下载
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp_146.hg38.vcf.gz & #snp库
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp_146.hg38.vcf.gz.tbi & #snp库
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz & #1000人类基因组计划数据
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.tbi & 
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.gz &#参考基因组 
nohup wget -c ftp://gsapubftp- anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.fai & #参考基因组索引
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.dict & 
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz & 
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi & 
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/funcotator/funcotator_dataSources.v1.6.20190124s.tar.gz & #功能注释
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/bwa_index/gatk_hg38* &#bwa索引
## 人类外显子坐标bed文件
wget ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_human/CCDS.current.txt

cat CCDS.current.txt |grep  "Public" |  perl -alne '{/[(.*?)]/;next unless $1;$gene=$F[2];$exons=$1;$exons=~s/s//g;$exons=~s/-/ /g;print "$F[0] $_  $gene" foreach split/,/,$exons;}' | sort -u | bedtools sort -i |awk '{print "chr"$0"    0   +"}'  > hg38.exon.bed ###perl不理解，后面再去查资料

因为时间关系，小郭老师已经下载好，我们直接调用就可以了。原始数据用的肖医师的数据，直接软连接到自己的目录

###软连接到自己的目录
cd ~/project/1.raw_fq
ln -s /tmp/N163/*.fq.gz ./&&ln -s /tmp/T163/*.fq.gz ./
ls -lh#查看一下是否成功

链接名称与原文名称不一致的原因是经过改名操作，步骤如下

###改名字
cat ls *.fq.gz >rename.txt
cat rename.txt
###结果如下
N163_DHE15981-A84-A67_AHNL7WCCXY_L1_1.clean.fq.gz   
N163_DHE15981-A84-A67_AHNL7WCCXY_L1_2.clean.fq.gz   
T163_DHE15976-A73-A2_AHNL7WCCXY_L1_1.clean.fq.gz    
T163_DHE15976-A73-A2_AHNL7WCCXY_L1_2.clean.fq.gz    
T163_DHE15976-A73-A2_AHNL7WCCXY_L2_1.clean.fq.gz    
T163_DHE15976-A73-A2_AHNL7WCCXY_L2_2.clean.fq.gz    
T163_DHE15976-A73-A2_AHNL7WCCXY_L3_1.clean.fq.gz    
T163_DHE15976-A73-A2_AHNL7WCCXY_L3_2.clean.fq.gz    
cat ls *.fq.gz ｜cut -d "_" -f 1,4-5|sed "s/clean.//g" >new.txt
cat new.txt
###结果如下
N163_L1_1.fq.gz
N163_L1_2.fq.gz
N163_L2_1.fq.gz
N163_L2_2.fq.gz
N163_L3_1.fq.gz
N163_L3_2.fq.gz
T163_L1_1.fq.gz
T163_L1_2.fq.gz
T163_L2_1.fq.gz
T163_L2_2.fq.gz
T163_L3_1.fq.gz
T163_L3_2.fq.gz
paste rename.txt new.txt >rename1.txt
cat rename1.txt
###结果如下
N163_DHE15981-A84-A67_AHNL7WCCXY_L1_1.clean.fq.gz   N163_L1_1.fq.gz
N163_DHE15981-A84-A67_AHNL7WCCXY_L1_2.clean.fq.gz   N163_L1_2.fq.gz
N163_DHE15981-A84-A67_AHNL7WCCXY_L2_1.clean.fq.gz   N163_L2_1.fq.gz
N163_DHE15981-A84-A67_AHNL7WCCXY_L2_2.clean.fq.gz   N163_L2_2.fq.gz
N163_DHE15981-A84-A67_AHNL7WCCXY_L3_1.clean.fq.gz   N163_L3_1.fq.gz
N163_DHE15981-A84-A67_AHNL7WCCXY_L3_2.clean.fq.gz   N163_L3_2.fq.gz
T163_DHE15976-A73-A2_AHNL7WCCXY_L1_1.clean.fq.gz    T163_L1_1.fq.gz
T163_DHE15976-A73-A2_AHNL7WCCXY_L1_2.clean.fq.gz    T163_L1_2.fq.gz
T163_DHE15976-A73-A2_AHNL7WCCXY_L2_1.clean.fq.gz    T163_L2_1.fq.gz
T163_DHE15976-A73-A2_AHNL7WCCXY_L2_2.clean.fq.gz    T163_L2_2.fq.gz
T163_DHE15976-A73-A2_AHNL7WCCXY_L3_1.clean.fq.gz    T163_L3_1.fq.gz
T163_DHE15976-A73-A2_AHNL7WCCXY_L3_2.clean.fq.gz    T163_L3_2.fq.gz
###替换名字
cat rename.txt |while read id
do
arr=(${id})
fq=${arr[0]}
new=${arr[1]}
mv ${fq} ${new}
done
ls *.fq.gz

建立一个contig，方便后面用

###建立contig
cd ..
ls *.gz|cut -d "_" -f 1-2|uniq > ../config
cat contig
N163_L1
T163_L1
T163_L2
T163_L3

建立bwa索引

###构建索引
bwa index   -a bwtsw  \ #用bwtsw算法
-p bwa_index   Homo_sapiens_assembly38.fasta  \ #索引前缀已经参考基因组
1>hg38.bwa_index.log 2>&1   & #成功与否都输入到hg38.bwa_index.log

同样也是节约空间时间，直接用小郭老师构建好的索引

wes流程

格式转换+质量控制

开始正式的流程,因为下载好的文件是fq.gz格式，所以可以跳过格式转换这一步

###格式转换加质量控制
cat >wget.sh
#!/bin/bash
### 从sra到gatk #跳过，因为这一次不需要
cat config |  while read id
do 

###sra2fq 
    # if [ ! -f ./1.raw_fq/${id}_1.fastq.gz ]
    # then
    # time fastq-dump --gzip --split-3 -O ./1.fq_raw 0.sra/${id}.sra
    # fi
    
### trim # 质量控制
    if [ ! -f ./3.clean_fq/${id}_1_val_1.fq.gz ]
    then
    fq1=./1.raw_fq/${id}_1.fq.gz
    fq2=./1.raw_fq/${id}_2.fq.gz
    time trim_galore -q 25 --phred33 --length 36 -e 0.1 --stringency 3 --gzip --paired -o ./3.clean_fq $fq1 $fq2 #参数含义：q 去掉q25以下的序列；--length去除长度小于36的序列；--phred 33 or 64；e 错误率 0.1； --stringency Overlap with adapter sequence required to trim a sequence，一般为3就行了； --gzip 输出是gz格式；--paired 双末端； -o 输出位置
    fi 
    done
cat wget.sh #检查脚本
nohup bash wget.sh &

*这里学到新的写脚本的方式，用if命令批量下载，如果不存在就下载，存在就跳过。

zcat N163_L1_1_val_1.fq.gz|wc -l
48485276
zcat N163_L1_1_val_2.fq.gz|wc -l
48485276

可以看到数目是一致的，原因是因为序列不同，所以压缩率有差异！

质量检测

一起调用fast_qc，速度更快

cat >wget1.sh
#!/bin/bash
cat config |  while read id
do 
###raw_qc 
    time fastqc -t 2 ./1.raw_fq/${id}*fq.gz -o ./2.qc/raw_qc
###clean_qc 
    time fastqc -t 2 ./3.clean_fq/${id}*fq.gz -o ./2.qc/clean_qc
done
cat wget1.sh #检查脚本
nohup bash wget1.sh &

质量控制结果对比

序列比对

选用bwa进行mapping

cat >wget2.sh
#!/bin/bash
INDEX=/home/hcguo/data/gatk/bwa_index
cat config |  while read id
do 
if [ ! -f ./4.align/${id}.bam ]
    then
    echo "start align for ${id}" `date`
    fq1=./3.clean_fq/${id}_1_val_1.fq.gz
    fq2=./3.clean_fq/${id}_2_val_2.fq.gz
    bwa mem -M -t 2 -R "@RG\tID:$id\tSM:$id\tLB:WXS\tPL:Illumina"  ${INDEX}  $fq1 $fq2 | samtools sort -@ 2 -m 1G  -o  ./4.align/${id}.bam -
#参数含义分别是：M，处理同一个reads比对到参考基因组上不同位置的情况，用来处理同一个reads比对到参考基因组上不同位置的情况；t 表示线程; R bam文件开头会有这一个header；索引；双末端；用samtools排序； t 线程； m 内存； o 输出位置  
    echo "end align for ${id}" `date`
    fi
    done
    cat wget2.sh #检查脚本
nohup   bash wget2.sh &
cd ~/project/4.align
ls -lh *bam|cut -d " " -f 5-
##结果如下
1.7G 12月 16 18:18 N163_L1.bam
3.5G 12月 16 20:55 T163_L1.bam

cat >wget3.sh
#!/bin/bash
cd 5.gatk
BAM=../4.align/${id}.bam
ref=/home/hcguo/data/gatk/Homo_sapiens_assembly38.fasta
bed=/home/hcguo/data//GRCh38.bed
GATK=/home/hcguo/biosoft/gatk-4.1.4.0/gatk
snp=/home/hcguo/data/gatk/dbsnp_146.hg38.vcf.gz
indel=/home/hcguo/data/gatk/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
cat config |  while read id
do 
if [ ! -f  ok.gatk.${id}.status ]
then
    echo "start gatk for ${id}" `date`
    $GATK --java-options "-Xmx4G -Djava.io.tmpdir=./" MarkDuplicates \ #去除重复
    -I ${BAM} \
    --REMOVE_DUPLICATES=true \ #去除重复
    -O ${id}_marked.bam \ #重矫正的bam
    -M ${id}.metrics \
    1>${id}_log.mark 2>&1 #节约空间
    
    samtools index -@ 4 -b ${id}_marked.bam

    $GATK --java-options "-Xmx4G -Djava.io.tmpdir=./"  BaseRecalibrator \ #碱基质量分数重校准
    -R ${ref}  \
    -I ${id}_marked.bam  \ #重矫正的bam
    --known-sites ${snp} \ #下载的snp库
    --known-sites ${indel} \ #下载的indel库
    -O ${id}_recal.table \ #输出为表格
    1>${id}_log.recal 2>&1 

    $GATK --java-options "-Xmx4G -Djava.io.tmpdir=./"   ApplyBQSR \ #应用碱基质量分数重校准
    -R ${ref}  \ #参考基因组
    -I ${id}_marked.bam  \ #重矫正的bam
    -bqsr ${id}_recal.table \ #上一步输出的表格
    -O ${id}_bqsr.bam \ #输出
    1>${id}_log.ApplyBQSR  2>&1 
##aplotypeCaller，简称HC，能过通过对活跃区域（也就是与参考基因组不同处较多的区域）局部重组装，同时寻找SNP和INDEL
#HC有一个GVCF模式，产生各个样本的中间基因组文件(gVCF),然后用于多样本的联合基因定型(joint genotyping)
    $GATK --java-options "-Xmx4G -Djava.io.tmpdir=./" HaplotypeCaller -ERC GVCF \
    -R ${ref} \
    -I ${id}_bqsr.bam \
    --dbsnp ${snp} \
    -L ${bed} \
    -O ${id}_raw.vcf \
    1>${id}_log.HC 2>&1 
    
    touch ok.gatk.${id}.status #创建空文件夹来计算程序运行数据
    echo "end gatk for ${id}" `date`
    sleep 100s

fi
cd ..
done
cat wget3.sh#检查脚本
nohup bash wget3.sh &

明天继续做剩下的部分包括snp与cnv已经结果的可视化！未完待续

下午 看医学方的视频 (作业)

看到群里的作业

想到几种不同的结果，

1.如果我想除去na的值

newdata<-na.omit(olddata)

2.如果我想将na的值换成一个特定的值，比如0

probes_expr[is.na(probes_expr)]<-0

3.换成每行中的均值，（根据上面的将0改成每一行的均值，但又的匹配每一行，所以谷歌了,没想到真有这样的包）

install.packages("Hmisc")
?imputed

probes_expr[is.na(probes_expr)<-(apply(probes_expr, 1, function(x){
  impute(x,mean)
}))

melt 函数数据拆分
dcast()函数合并

x<-data.frame(id=1:6,
              name=c("wang","zhang","li","chen","zhao","song"),
              shuxue=c(89,85,68,79,NA,53),
              yuwen=c(77,68,NA,87,92,63))
library(reshape2)
x1<-melt(x,id=c("id","name"))

结果如下格式


image.png

image.png

dcast输出结果为数据框。

dcast(x1,id~variable,mean,margins=T)

可以看到，边缘多了两列汇总数据是对行列求平均的结果


image.png

函数还是要多理解

©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 204,053评论 6赞 478
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 85,527评论 2赞 381
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 150,779评论 0赞 337
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,685评论 1赞 276
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,699评论 5赞 366
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,609评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 37,989评论 3赞 396
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,654评论 0赞 258
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 40,890评论 1赞 298
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,634评论 2赞 321
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,716评论 1赞 330
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,394评论 4赞 319
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 38,976评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,950评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,191评论 1赞 260
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 44,849评论 2赞 349
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,458评论 2赞 342

推荐阅读更多精彩内容

  • 大前端的三层思考-2019-12-21
    是否选择大前端？ 虽然现在大前端很热门，是不是适合自己呢？这个需要想一想 技术的广度和深度，哪个优先？ 可以简单类...
    勇往直前888阅读 347评论 0赞 3
  • java的数据类型和java命名规范----------2019-12-21
    今天继续昨天话题，谈谈java的数据类型和java命名规范 一。数据类型 java数据类型共有三类八种，具体如下图...
    lucky珂阅读 180评论 0赞 0
  • 2019-12-21
    1.整理校宝，统计名辉豪庭校区报名花名册数据并分组非陪给校区老师，安排续费任务（2.5H） 2.联系昨晚对接转介绍...
    科比黑曼巴阅读 115评论 0赞 0
  • 2019-07-13
    日精进打卡第470天 知～学 《六项精进 》0遍共399遍 《大学》0遍共390遍 《六项精进通篇》0遍共1 【经...
    常庚_c041阅读 164评论 0赞 0
  • 别谈成功 先照顾好你自己
    最近，由于季节变化，很多人都因为温度的变差而生病了。当然我这个多病体质也没有幸免逃脱，也华丽丽的生起病了。说实话，...
    Jessica夏沫阅读 281评论 0赞 0