非模式物种真的到处都是bug, 每次出bug我都觉得自己学艺不精，当然也的的确确是学艺不精。记录一下RNA速率分析踩过的各种大坑。
RNA速率分析需要输入的文件
· CellRanger 处理文件
· 基因组的注释文件
· 重复序列的注释文件
①CellRanger处理文件结果比较好得到，主要是Bam文件和barcodes.tsv。
②注释文件我本来以为很简单，都是在cellranger时候准备好的，但是真的大错特错。这部分后续再说。
③重复序列注释文件，一般组装基因组的时候会涉及到。但是也可以自己生成。
重复序列注释文件
重复序列注释方法（RepeatMasker）分为同源比对和从头预测两类。
①基于同源性的预测，用针对 RepBase（www.repeatmasker.org）和Dfam(https://www.dfam.org/home
)筛选硬骨鱼基因组中已知的转座因子。（但是这个库之前是免费的，现在好像已经不免费了，所以我用的是网上的2018版的）
②从头预测由 RepeatModeler 和 LTR？构建的重复文库，并据此库重新预测重复序列。
参考链接：https://www.jianshu.com/p/a0dee7c5fdef
https://www.biochen.org/cn/blog/2021/%E4%BD%BF%E7%94%A8RepeatModeler%E4%BB%8E%E5%A4%B4%E9%A2%84%E6%B5%8B%E5%9F%BA%E5%9B%A0%E7%BB%84%E9%87%8D%E5%A4%8D%E5%BA%8F%E5%88%97/
安装软件:这里需要非常多的依赖软件，主要是RepeatMasker和RepeatModeler这两个软件及其依赖软件。
（一）整合数据库
1. 进入RepeatMasker目录
wget https://github.com/chenyangkang/Repbase-Dfam/blob/main/RepBaseRepeatMaskerEdition-20181026.tar.gz
tar -zxvf RepBaseRepeatMaskerEdition-20181026.tar.gz
RepeatMasker/Libraries中应当包含两个文件RMRBSeqs.embl和README.RMRBSeqs
2.Dfam下载及合并
wget -c https://www.dfam.org/releases/Dfam_3.4/families/Dfam_curatedonly.embl.gz
wget -c https://www.dfam.org/releases/Dfam_3.4/families/Dfam_curatedonly.h5.gz
wget -c https://www.dfam.org/releases/Dfam_3.4/families/Dfam_curatedonly.hmm.gz
Curated（也就是不包含从头预测的物种）,因为repeatmodler可以进行从头预测，所以只下载Curated即可。
gunzip Dfam_curatedonly.embl.gz
gunzip Dfam_curatedonly.h5.gz
gunzip Dfam_curatedonly.hmm.gz
mv Dfam_curatedonly.embl Dfam.embl
mv Dfam_curatedonly.h5 Dfam.h5
mv Dfam_curatedonly.hmm Dfam.hmm
3.重新congifure RepeatMasker
./configure
combining Dfam+RepBase生成以及配置结束以后Dfam和RepBase的版本号，同时Libraries里会生成一个RepeatMaskerLib.h5(两个数据库的整合版本)
4.导出已有数据库中的近缘物种
famdb.py -i RepeatMaskerLib.h5 families \
-f embl -a -d Teleostei> teleost.embl #这里主要是根据然后根据taxonomy分类命名找
我主要找的是硬骨鱼纲，为了和RepeatModeler合并，要将文件转换为fasta格式。
buildRMLibFromEMBL.pl teleost.embl> teleost.fasta
buildRMLibFromEMBL.pl脚本在RepeatMasker/util文件夹里
（二）从头预测
###用序列文件构建数据库
BuildDatabase -name ${i} -engine rmblast ${i}.fna
###RepeatModeler根据数据库比对预测repeats
RepeatModeler -database ${i} -pa 5
其中${i}为基因组文件名
nohup RepeatModeler -database ${i} -pa 12 -LTRStruct > run.out &
顺利的话最后得到一个文件夹，以及一个${i}-families.fa，一个${i}-families.stk文件。前者就是从头预测得到的数据库
（三）两个数据库合并
数据库和从头预测的结果合并
cat ${i}-families.fa teleost.fasta > all_${i}_final.fasta
（四）RepeatMasker注释
RepeatMasker -xsmall -poly \
-pa 5 -lib all_${i}_final.fasta -engine ncbi ${i}.fna
比对结果
我尝试了三种方法，分别是数据库、从头预测、以及数据库和从头预测整合分析。三个的方法预测结果差异比较大，从头预测和整合分析结果差异比较小。但是我分析的是草鱼，我预测的重复序列50%多，我师兄说草鱼和团头鲂比较接近，他文章里面给的结果是40%多。但是具体是哪里不一样我也不太清楚。但就按照这个继续往后做。
image.png
本来以为重复序列做好就可以直接RNA velocyto了，但是，没想到前面有跟大的坑
RNA 速率分析
安装和配置velocyto，具体参见网上一些教程，这里想记录一下bug
image.png
跑了不到一分钟，然后exit了，打开error.log，好家伙，velocyto本质上是python运行的软件，出来直接是python的报错，我python学艺不精，直接把我整懵了。但是仔细一看推测注释文件的问题。注释文件能有什么问题？我去github上找，但是我看velocyto软件的文档里面写的注释文件遵循GENCODE的格式，但是我真的都差拿放大镜找不同了，但是也没有什么不同。
但功夫不负有心人，我找到了github上非模式物种关于RNA速率分析的一个问题
https://github.com/velocyto-team/velocyto.py/issues/73
其实主要就是基因注释的不够规范，最好能先用IGV检查基因组的注释。
我的问题主要是外显子区域发生重叠导致内含子小于零出现了报错。
发生这种情况最直接的办法就是保留该基因的最长的外显子，删掉其余外显子。删除重叠外显子的脚本
#首先对注释文件进行排序
cat genome_mt.gtf|awk '$3=="exon"{print $0}'|sort -t $'\t' -k 1,1 -k 7,7 -k 4n,4 -k 5n,5 >sort.gtf
#剔除重叠的短的外显子
cat genome_mt.gtf|awk '$3=="exon"{print $0}'|sort -t $'\t' -k 1,1 -k 7,7 -k 4n,4 -k 5n,5 >sort.gtf
awk 'BEGIN{
RS="\n"
FS="\t"
flag=1
}
$3=="exon"{
if(!start[$1","$7]){
start[$1","$7]=$4
end[$1","$7]=$5
flag=1
print $0
}
else{
s=start[$1","$7]
e=end[$1","$7]
if((s<=$4&&$4<=e)||(s<=$5&&$5<=$e)){
flag=0
}
else{
start[$1","$7]=$4
end[$1","$7]=$5
flag=1
print $0
}
}
}
$3!="exon"&&flag==1{print $0}
' sort.gtf > filtersort.gtf
这里又恶补了awk，我真的沙子堆的基础·······
然后进入CellRanger处理的文件中，进行velocyto
cellranger_outDir=~/sc/Grasscarp/BL/Day16_Control_BL
cellranger_gtf=~/genome/Grasscarp/database/filtersort.gtf
rmsk_gtf=~/genome/Grasscarp/repeatmask/denovo/genome_mt.fasta.out.gff
ls -lh $rmsk_gtf $cellranger_outDir $cellranger_gtf
velocyto run10x -m $rmsk_gtf $cellranger_outDir $cellranger_gtf
成功 🧡🧡🧡🧡🧡🧡🧡，撒花~~~~~~~~
image.png
网友评论