参考文献
[PMID: 27694136]
DRETools - F1000Research
数据下载
使用本实验室一个RNA-seq（GSE56087）和一个circRNA-seq（未发布）的数据作为discovery dataset 和 validation dataset。
RNA-seq数据下载
shell命令批量下载
#!/bin/sh
#调用prefetch下载SRA数据库的测序文件，并用pfastq-dump转化为fastq文件。
#需要在当前目录的download_list.txt文件中把SRA编号，每个一行事先准备好，具体格式参考example.txt
#建议使用NCBI的SRA RUN SELECTOR直接导出。
read -p "请输入下载列表文件名:" downloadlist
read -p "请选择单端(SE)，双端(PE)测序:" class
read -p "请输入要使用的线程数:" thread
if [ $class = "SE" ]; then
cat $downloadlist | while read line
do
echo $line
prefetch $line -o ./${line}.sra
pfastq-dump --threads $thread --outdir ./ ./${line}.sra
done
elif [ $class = "PE" ]; then
cat $downloadlist | while read line
do
echo $line
prefetch $line -o ./${line}.sra
pfastq-dump --split-3 --threads $thread --outdir ./ ./${line}.sra
done
else echo "错误的输入！"
fi
相应的accession number
SRR1200850
SRR1200851
SRR1200852
SRR1200853
SRR1200854
SRR1200855
SRR1200856
SRR1200857
SRR1200860
SRR1200861
SRR1200858
SRR1200859
SRR1200862
SRR1200863
SRR1200864
SRR1200865
SRR1200866
SRR1200867
SRR1200868
SRR1200869
SRR1200870
SRR1200871
SRR1200872
SRR1200873
SRR1200874
SRR1200875
SRR1200876
SRR1200877
SRR1200878
SRR1200879
SRR1200880
SRR1200881
SRR1200882
SRR1200883
SRR1200884
SRR1200885
数据分析
主要使用RNAEditor软件分析，环境配置按照官网的配置，注意相应软件的版本，命令需要在RNAEditor目录下运行，并且要求不能在SSH界面下运行，会报错（无法画图），需要在图形界面下打开终端运行。
shell批量分析
for ((i=1200851;i<1200886;i++))
do
fastp -w 16 -i ~/RNASEQ/EC_lina/SRR${i}_1.fastq -I ~/RNASEQ/EC_lina/SRR${i}_2.fastq -o ~/RNASEQ/EC_lina/SRR${i}_1_fastpedited.fastq -O ~/RNASEQ/EC_lina/SRR${i}_2_fastpedited.fastq
python RNAEditor.py -i ~/RNASEQ/EC_lina/SRR${i}_1_fastpedited.fastq ~/RNASEQ/EC_lina/SRR${i}_2_fastpedited.fastq -c configuration.txt
done
配置文件configuration.txt如下
#This file is used to configure the behaviour of RNAeditor
#Standard input files
refGenome = /home/zhou/rnaEditor_annotations/human/GRCH38/Homo_sapiens.GRCh38.dna.primary_assembly.fa
gtfFile = /home/zhou/rnaEditor_annotations/human/GRCH38/Homo_sapiens.GRCh38.83.gtf
dbSNP = /home/zhou/rnaEditor_annotations/human/GRCH38/dbSNP.vcf
hapmap = /home/zhou/rnaEditor_annotations/human/GRCH38/HAPMAP.vcf
omni = /home/zhou/rnaEditor_annotations/human/GRCH38/1000GenomeProject.vcf
esp = /home/zhou/rnaEditor_annotations/human/GRCH38/ESP_filtered
aluRegions = /home/zhou/rnaEditor_annotations/human/GRCH38/repeats.bed
output = default
sourceDir = /usr/local/bin/
maxDiff = 0.04
seedDiff = 2
standCall = 0
standEmit = 0
edgeDistance = 3
paired = True
keepTemp = True
overwrite = False
threads = 23
开23个核，每个样本分析时间大概在十个小时左右。
结果解析
结果实例
image.png
VCF文件
VCF文件包含所有的编辑站点以供进一步分析.
image.png
GCF文件
GCF文件保存了每个编辑站点的附加信息，如基因名称、片段、总读数、编辑读数和编辑比。
image.png
summary文件
summary文件显示每个基因的RNA编辑所处基因位置（如3‘UTR，5‘UTR）的数量。
image.png
bed文件
Editing island的bed文件，可供后续可视化等分析。
image.png
打包除了bam，sam，sai这两种文件以外的所有结果
tar -cvf RNAEditor_1.tar ~/RNASEQ/EC_lina/rnaEditor/ --exclude=*bam --exclude=*sam --exclude=*sai
pigz RNAEditor_1.tar
circRNA-seq数据集的验证工作
shell批量分析
#!/bin/sh
#把fastq文件路径存入数组
c=0
for file in `find ~/RNASEQ/CircularRNA-seq/rawdata/ -name *fastq.gz -print`
do
filelist[$c]=$file
((c++))
done
for ((i=0;i<${#filelist[@]};i=i+2))
do
tmp1=$(echo ${filelist[$i]}|sed 's/\.fastq\.gz/\_fastpedited\.fastq/g')
tmp2=$(echo ${filelist[$i+1]}|sed 's/\.fastq\.gz/\_fastpedited\.fastq/g')
#fastp预处理
fastp -w 16 -i ${filelist[$i]} -I ${filelist[$i+1]} -o $tmp1 -O $tmp2
#pigz -d $tmp1
#pigz -d $tmp2
python RNAEditor.py -i $tmp1 $tmp2 -c configuration.txt
done
configuration.txt不变
打包除了bam，sam，sai这两种文件以外的所有结果
tar -cvf RNAEditor_2.tar ~/RNASEQ/CircularRNA-seq/rawdata/Sample_ZYF-*/rnaEditor --exclude=*bam --exclude=*sam --exclude=*sai
pigz RNAEditor_2.tar
后续分析
DREtools差异分析
#Merge editing sites from multiple samples
dretools edsite-merge --min-editing 3 --min-coverage 5 --min-samples 3 --vcf ./*/*.editingSites.vcf > DRE/consensus_sites.vcf
#Calculate sample-wise EPK
for i in `find ./SRR*/ -name *fastpedited.bam -print` ;
do
name=$(echo ${i}|sed 's/_1.*bam//'|sed 's/\.\///');
dretools sample-epk --name $name --vcf ./DRE/consensus_sites.vcf --alignment $i > ./DRE/${name}.sample_epk.tsv;
done
#Calculate site-wise EPK
for i in `find ./SRR*/ -name *fastpedited.bam -print` ;
do
name=$(echo ${i}|sed 's/_1.*bam//'|sed 's/\.\///');
dretools edsite-epk --vcf ./DRE/consensus_sites.vcf --alignment $i > ./DRE/${name}.edsite_epk.tsv;
done
#Find differentially edited editing sites
dretools edsite-diff --max-depth-cov 5.0 --min-depth 2 \
--names normal,tumor \
--sites ./DRE/consensus_sites.vcf \
--sample-epk ./DRE/SRR1200850.sample_epk.tsv,./DRE/SRR1200851.sample_epk.tsv,./DRE/SRR1200854.sample_epk.tsv,./DRE/SRR1200855.sample_epk.tsv,./DRE/SRR1200858.sample_epk.tsv,./DRE/SRR1200859.sample_epk.tsv,./DRE/SRR1200862.sample_epk.tsv,./DRE/SRR1200863.sample_epk.tsv,./DRE/SRR1200866.sample_epk.tsv,./DRE/SRR1200867.sample_epk.tsv,./DRE/SRR1200872.sample_epk.tsv,./DRE/SRR1200873.sample_epk.tsv,./DRE/SRR1200874.sample_epk.tsv,./DRE/SRR1200875.sample_epk.tsv,./DRE/SRR1200878.sample_epk.tsv,./DRE/SRR1200879.sample_epk.tsv,./DRE/SRR1200882.sample_epk.tsv,./DRE/SRR1200883.sample_epk.tsv \
./DRE/SRR1200852.sample_epk.tsv,./DRE/SRR1200853.sample_epk.tsv,./DRE/SRR1200856.sample_epk.tsv,./DRE/SRR1200857.sample_epk.tsv,./DRE/SRR1200860.sample_epk.tsv,./DRE/SRR1200861.sample_epk.tsv,./DRE/SRR1200864.sample_epk.tsv,./DRE/SRR1200865.sample_epk.tsv,./DRE/SRR1200868.sample_epk.tsv,./DRE/SRR1200869.sample_epk.tsv,./DRE/SRR1200870.sample_epk.tsv,./DRE/SRR1200871.sample_epk.tsv,./DRE/SRR1200876.sample_epk.tsv,./DRE/SRR1200877.sample_epk.tsv,./DRE/SRR1200880.sample_epk.tsv,./DRE/SRR1200881.sample_epk.tsv,./DRE/SRR1200884.sample_epk.tsv,./DRE/SRR1200885.sample_epk.tsv \
--site-epk ./DRE/SRR1200850.edsite_epk.tsv,./DRE/SRR1200851.edsite_epk.tsv,./DRE/SRR1200854.edsite_epk.tsv,./DRE/SRR1200855.edsite_epk.tsv,./DRE/SRR1200858.edsite_epk.tsv,./DRE/SRR1200859.edsite_epk.tsv,./DRE/SRR1200862.edsite_epk.tsv,./DRE/SRR1200863.edsite_epk.tsv,./DRE/SRR1200866.edsite_epk.tsv,./DRE/SRR1200867.edsite_epk.tsv,./DRE/SRR1200872.edsite_epk.tsv,./DRE/SRR1200873.edsite_epk.tsv,./DRE/SRR1200874.edsite_epk.tsv,./DRE/SRR1200875.edsite_epk.tsv,./DRE/SRR1200878.edsite_epk.tsv,./DRE/SRR1200879.edsite_epk.tsv,./DRE/SRR1200882.edsite_epk.tsv,./DRE/SRR1200883.edsite_epk.tsv \
./DRE/SRR1200852.edsite_epk.tsv,./DRE/SRR1200853.edsite_epk.tsv,./DRE/SRR1200856.edsite_epk.tsv,./DRE/SRR1200857.edsite_epk.tsv,./DRE/SRR1200860.edsite_epk.tsv,./DRE/SRR1200861.edsite_epk.tsv,./DRE/SRR1200864.edsite_epk.tsv,./DRE/SRR1200865.edsite_epk.tsv,./DRE/SRR1200868.edsite_epk.tsv,./DRE/SRR1200869.edsite_epk.tsv,./DRE/SRR1200870.edsite_epk.tsv,./DRE/SRR1200871.edsite_epk.tsv,./DRE/SRR1200876.edsite_epk.tsv,./DRE/SRR1200877.edsite_epk.tsv,./DRE/SRR1200880.edsite_epk.tsv,./DRE/SRR1200881.edsite_epk.tsv,./DRE/SRR1200884.edsite_epk.tsv,./DRE/SRR1200885.edsite_epk.tsv \
> ./DRE/diff_sites.tsv
#Detect editing islands
dretools find-islands \
--min-editing 3 \
--min-coverage 5 \
--min-length 20 \
--min-points 5 \
--epsilon 50 \
--vcf ./*/*.editingSites.vcf > ./DRE/islands.bed
#Calculate Island-EPK
for i in `find ./SRR*/ -name *fastpedited.bam -print` ;
do
name=$(echo ${i}|sed 's/_1.*bam//'|sed 's/\.\///');
dretools region-epk --regions ./DRE/islands.bed --vcf ./DRE/consensus_sites.vcf --alignment $i > ./DRE/${name}.island_epk.tsv;
done
#Find differentially edited islands
dretools region-diff --regions ./DRE/islands.bed --min-area 1 --min-depth 1 \
--names normal,tumor \
--sample-epk ./DRE/SRR1200850.sample_epk.tsv,./DRE/SRR1200851.sample_epk.tsv,./DRE/SRR1200854.sample_epk.tsv,./DRE/SRR1200855.sample_epk.tsv,./DRE/SRR1200858.sample_epk.tsv,./DRE/SRR1200859.sample_epk.tsv,./DRE/SRR1200862.sample_epk.tsv,./DRE/SRR1200863.sample_epk.tsv,./DRE/SRR1200866.sample_epk.tsv,./DRE/SRR1200867.sample_epk.tsv,./DRE/SRR1200872.sample_epk.tsv,./DRE/SRR1200873.sample_epk.tsv,./DRE/SRR1200874.sample_epk.tsv,./DRE/SRR1200875.sample_epk.tsv,./DRE/SRR1200878.sample_epk.tsv,./DRE/SRR1200879.sample_epk.tsv,./DRE/SRR1200882.sample_epk.tsv,./DRE/SRR1200883.sample_epk.tsv \
./DRE/SRR1200852.sample_epk.tsv,./DRE/SRR1200853.sample_epk.tsv,./DRE/SRR1200856.sample_epk.tsv,./DRE/SRR1200857.sample_epk.tsv,./DRE/SRR1200860.sample_epk.tsv,./DRE/SRR1200861.sample_epk.tsv,./DRE/SRR1200864.sample_epk.tsv,./DRE/SRR1200865.sample_epk.tsv,./DRE/SRR1200868.sample_epk.tsv,./DRE/SRR1200869.sample_epk.tsv,./DRE/SRR1200870.sample_epk.tsv,./DRE/SRR1200871.sample_epk.tsv,./DRE/SRR1200876.sample_epk.tsv,./DRE/SRR1200877.sample_epk.tsv,./DRE/SRR1200880.sample_epk.tsv,./DRE/SRR1200881.sample_epk.tsv,./DRE/SRR1200884.sample_epk.tsv,./DRE/SRR1200885.sample_epk.tsv \
--region-epk ./DRE/SRR1200850.island_epk.tsv,./DRE/SRR1200851.island_epk.tsv,./DRE/SRR1200854.island_epk.tsv,./DRE/SRR1200855.island_epk.tsv,./DRE/SRR1200858.island_epk.tsv,./DRE/SRR1200859.island_epk.tsv,./DRE/SRR1200862.island_epk.tsv,./DRE/SRR1200863.island_epk.tsv,./DRE/SRR1200866.island_epk.tsv,./DRE/SRR1200867.island_epk.tsv,./DRE/SRR1200872.island_epk.tsv,./DRE/SRR1200873.island_epk.tsv,./DRE/SRR1200874.island_epk.tsv,./DRE/SRR1200875.island_epk.tsv,./DRE/SRR1200878.island_epk.tsv,./DRE/SRR1200879.island_epk.tsv,./DRE/SRR1200882.island_epk.tsv,./DRE/SRR1200883.island_epk.tsv \
./DRE/SRR1200852.island_epk.tsv,./DRE/SRR1200853.island_epk.tsv,./DRE/SRR1200856.island_epk.tsv,./DRE/SRR1200857.island_epk.tsv,./DRE/SRR1200860.island_epk.tsv,./DRE/SRR1200861.island_epk.tsv,./DRE/SRR1200864.island_epk.tsv,./DRE/SRR1200865.island_epk.tsv,./DRE/SRR1200868.island_epk.tsv,./DRE/SRR1200869.island_epk.tsv,./DRE/SRR1200870.island_epk.tsv,./DRE/SRR1200871.island_epk.tsv,./DRE/SRR1200876.island_epk.tsv,./DRE/SRR1200877.island_epk.tsv,./DRE/SRR1200880.island_epk.tsv,./DRE/SRR1200881.island_epk.tsv,./DRE/SRR1200884.island_epk.tsv,./DRE/SRR1200885.island_epk.tsv \
> ./DRE/diff_islands.tsv
DREtools结果展示
image.png
MET-DB交集
