用picard制作外显子坐标文件

---------------

Nickier

2018-12-07

---------------

第一次用picard，不知道怎么用，先用conda安装了一个

# 安装前先search一下
conda search picard

# 发现有很多个版本
# Name                 Version           Build Channel            
# picard                     1.56               0 anaconda/cloud/bioconda
# picard                     1.56               1 anaconda/cloud/bioconda
# picard                     1.97               0 anaconda/cloud/bioconda

# 那就随便下一个吧
conda install -y picard

# 然后调出帮助文档
picard --help

找到需要的命令BedToIntervalList

BedToIntervalList.jpeg

然后用picard BedToIntervalList看看有什么参数，发现

# Usage example:
#
# java -jar picard.jar BedToIntervalList \
# I=input.bed \
# O=list.interval_list \
# SD=reference_sequence.dict

而我的input.bed是hg38.chr.bed，reference_sequence.dict是/public/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.dict

如果没有bed文件，就到CCDS数据库下载

# 下载方法
wget ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_human/CCDS_exons.20180614.txt

# 然后用perl转换格式
cat CCDS_exons.20180614.txt |perl -alne  '{/\[(.*?)\]/;next unless $1;$gene=$F[2];$exons=$1;$exons=~s/\s//g;$exons=~s/-/\t/g;print "$F[0]\t$_\t$gene" foreach split/,/,$exons;}'|sort -u |bedtools sort  -i >exon_probe.hg38.gene.bed

# 再提取
cat exon_probe.hg38.gene.bed|awk '{print "chr"$0}' >hg38.chr.bed

​
有了bed文件之后，就可以用一下命令啦

java -jar picard.jar BedToIntervalList \

I=hg38.chr.bed \

O=list.interval_list \

SD=/public/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.dict

然而最不想发生的事情发生，报错：Error: Unable to access jarfile picard.jar

544187879816.png

查了一下解决报错的方法，原来java需要调用.jar的全路径。

java -jar /home/hcguo/miniconda3/envs/wes/share/picard-2.18.17-0/picard.jar BedToIntervalList \

I=hg38.chr.bed \

O=list.interval_list \

SD=/public/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.dict

这样就解决问题啦，也达到需求了。生成了list.interval_list文件

## /public/biosoft/GATK/resources/bundle/hg38
## bed to intervals_list
## /home/hcguo/HT_1/tmp/hg38.chr.bed

GATK=/home/jmzeng/biosoft/gatk4/gatk-4.0.6.0/gatk
java -jar picard.jar BedToIntervalList \
I=hg38.chr.bed \
O=list.interval_list \
SD=/public/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.dict

## Preprocess Intervals

$GATK PreprocessIntervals \
-L list.interval_list \
--sequence-dictionary /public/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.dict \
--reference /public/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta \
--padding  250 \
--bin-length  0 \
--interval-merging-rule OVERLAPPING_ONLY \
--output targets.preprocessed.interval.list