在「Bionano系列」光学图谱混合组装应该怎么做?这篇文章中,我展示了下面这张图。 和之前的图不同的是,我加了几个箭头,这些箭头所指向的区域的特征就是,这些区域并未被Bionano所覆盖。如果不去思考这些区域到底是什么,直接进行混合组装,那么这其实对最后结果的不负责任。因为这完全可能是组装软件没有正确的处理错误的overlap,将不应该连接的序列连接在一起(尽管这个概率不高)。 我的直观猜测就是,这些区域应该是重复序列区域。毕竟Bionano标记技术依赖于酶识别特定位点进行酶切加上荧光标记,重复序列要么会因酶切密度太高,相机的分辨率达不到,而识别失败,要么是酶切位点过少,信号太弱。 那么我应该如何验证这个猜想?通过几天的文献翻阅和尝试,我用重复序列数量和基因数量的相对比值进行衡量。 命令行的代码如下(没有考虑文件的相对位置) # 利用拟南芥的原本CDS进行注释 gmap_build -D index -d R05C0144 ../R05C0144.fa & gmap -t 20 -D index -d R05C0144 -f gff3_gene ../Athaliana_cds.fa > cds_gene.gff3 2> log.txt & # 重复序列注释 RepeatMasker -e ncbi -species arabidopsis -pa 30 -gff -dir . ../R05C0144.fa & # GFF转成BED awk 'BEGIN{OFS="\t"} {print $1,$4,$5}' ../repeat_annotation/R05C0144.fa.out.gff > repeat.bed grep -w 'gene' ../gene_annotation/cds_gene.gff3| awk 'BEGIN{OFS="\t"} {print $1,$4,$5}' | bedtools sort -i - > gene.bed bedtools makewindows -w 100000 -g ../R05C0144.txt > windows_100k.bed bedtools coverage -a windows_100k.bed -b repeat.bed > repeat_stat.bed bedtools coverage -a windows_100k.bed -b gene.bed > gene_stat.bed R代码如下 gene_df <- read.table("R05C0144/feature_stat/gene_stat.bed", sep="\t", stringsAsFactors = F) repeat_df <- read.table("R05C0144/feature_stat/repeat_stat.bed", sep="\t", stringsAsFactors = F) options(scipen=999) contig <- "contig2" repeat_ctg <- repeat_df[repeat_df$V1 == contig,] gene_ctg <- gene_df[gene_df$V1 == contig,] combine_df <- data.frame(pos=(repeat_ctg$V2 + repeat_ctg$V3) / 2, repeat_num=repeat_ctg$V4, gene_num=gene_ctg$V4) combine_df$total = combine_df$repeat_num + combine_df$gene_num combine_df$gene_ratio <- combine_df$gene_num / combine_df$total * 100 combine_df$repeat_ratio <- combine_df$repeat_num / combine_df$total * 100 plot(combine_df$pos, combine_df$gene_ratio, type="l", ylim=c(0,100), xlab="position", ylab="percent", col="blue") lines(combine_df$pos, combine_df$repeat_ratio, col="red") abline(v=7.85*1e6) 我检查了一些区间,的确是重复序列比例高于基因比例,当然还有一些区间不是。说明重复序列并不是光学图谱未覆盖的主要原因。 当然对于拟南芥这种有着高质量基因组的物种而言,我们还可以进行共线性分析。不过对于这些N50在4M左右,而且低杂合的基因组,其实都不需要太操心这种错误。 我这里也就验证了一种可能性,后续还得检查了一下其他原因,说不定仅仅是光学图谱的深度不够而已。

AutoNoise + SplitBNX: 这一步会将bnx和参考的cmap文件进行比对,估算出噪声系数,然后把bnx进行拆分便与后续比对 Pairwse: 这一步进行molecules之间的两两比较,寻找overlap, 结果存放在"align"文件夹下 Assembly: 根据两两比对结果,通过OLC算法进行组装,结果在"contigs/exp_unrefined"下,合并后的文件为"EXP_UNREFINED.cmap", 同时还会将"EXP_UNREFINED.cmap"和参考基因组的cmap进行比对,结果放在"contigs/exp_unrefined/exp_unrefined/alignref", 此外还将拆分后Bnx文件和参考基因组的cmap文件进行比对,结果放在"contigs/alignmolvref"下 refineA: 第3步得到的图谱先会进行第一次优化 输出结果在"contigs/exp_refineA", 这一步不会使用所有的原始数据,而是pairwise阶段用的质量比较好的分子,所以速度会快一些 refineB: 在第4步的基础上进行第二次的优化, 输出结果在"contigs/exp_refineB0"和"contigs/exp_refineB1". 这一步会使用所有的输入原始数据,速度稍微慢一些。第一轮和第二轮的结果会将地覆盖度的区域进行打断,然后更新标记和标记的位置。 Extension and merge: 将上一步的contig回贴到参考基因组的map,进行延伸和合并,这一步可以迭代3-5次。中间结果在"contigs/exp_extensionX_X"和"contigs/exp_mrgX" Final refinement: 最后一步的优化 SV Decetion: 在有参考基因组的前提下,最后还会寻找一些大规模的结构变异。 上述这些流程由Solve工具中Pipeline文件夹下的脚本pipelineCL.py进行控制,以之前的过滤后的人类数据为例 python /opt/biosoft/Solve3.3_10252018/Pipeline/10252018/pipelineCL.py \ -T 96 -N 4 -f 1 \ -i 5 \ # 延伸和合并的迭代次数,默认是1,介于0~20之间 -b molecules120k.bnx \ # 输入的BNX文件 -r NA12878_CTTAAG_0kb_0labels.cmap \ #参考的reference,可选 -l Assembly \ # 输出文件夹 -y \ # 自动确定噪声参数,需要-r -t /opt/biosoft/Solve3.3_10252018/RefAligner/7915.7989rel \ # RefAligner和Assembler的所在文件夹 -a /opt/biosoft/Solve3.3_10252018/RefAligner/7915.7989rel/optArguments_nonhaplotype_saphyr_human.xml # 配置参数文件 上面的-a参数最为重要,因为指定的xml文件控制了流程中每一步的具体参数,所以要慎重选择。 重点1: XML文件命名解释说明: irys/saphyr: 数据来源仪器 DLE1: DEL1标记系统 human: 物种是人类 BG: big genome. 大于5G,主要是优化内存使用 noES: no extend and split: 即便等位基因里有超过30kbp的结构变异,也不要将他们分开 haplotype/nonhaplotype: 单倍型优化指的是将那些含有超过500bp或者更大的SV差异的等位基因进行分开,不推荐用于非人类基因组组装是用haplotype,这会导致组装结果过度碎片化,组装的基因组会变大. nocut: 不对CMPR(complex multipath regions)进行拆分,所谓的CMPR指的是长度超过130kbp的高度重复序列,因为相似度过高,组装的时候不知道如何处理。 对于非人类的物种,推荐参数为: 除非是小鼠的SV缺失,大部分情况都用nonhaplotype,用于后续的Hybrid scaffold Pipeline(HS) 要noES 是否cut看情况而定。大部分人喜欢不cut。 对于人类, 推荐参数为: 仅在SV检测时用haplotype, 对于HS用nonhaplotype 大部分情况下加上CMPR, 除非你知道在CMPR区间上有SV,才使用nocut 重点2: 如何设置组装时pairwise alignment中的-T参数。 基本原则是: 标记每增加一个,p值降低100倍; 基因组每增加一个数量级,p值降低100倍。 对于大于1Gb, 标记密度小于 15/100bkp的情况,设置为1e-11, 对于基因组大于1Gb, 标记密度大于15/100 kbp,每增加一个标记,就降低100倍,例如17/100kbp, 推荐1e-15 此外可以用-B跳过部分流程,-e则是输出文件的前缀,-x表示在自动去噪后退出。 如果有集群可用参数-C。如果基因组比较差, 可用-R 参数进行初步组装,然后基于第一个版本进一步的组装。还有一些参数用于关闭一些默认启动的功能 -A: 不将bnx文件和最后的contig比较 -m: 不将bnx和参考基因组比较 -E: 不检查输出的完整性 运行过程中, 可以用"grep 'Executing' bionanoAssembly/exp_pipelineReport.txt" 查看执行进度 最后输出结果是"contigs/exp_refineFinal1/EXP_REFINEFINAL1.cmap",而结果好坏则要看"exp_informaticsReportSimple.txt", 两个核心标准 Bionano图谱应该占原来的物理图谱的90%以上 Bionano图谱的N50 会由于物种不同有很大差异,动物一般在1M以上,植物不确定。DLE系统差异更明显

从这部分开始,就开始涉及一些软件的操作和数据分析,因此在进入正文之前,我们需要准备好环境。 第一步:从 https://bionanogenomics.com/library/datasets/下载人类测试数据集,以及对应的NA12878人类基因组。 wget http://bnxinstall.com/publicdatasets/DLS/20180413_NA12878_S3_compressed.tar.gz tar xf 20180413_NA12878_S3_compressed.tar.gz wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/002/077/035/GCA_002077035.3_NA12878_prelim_3.0/GCA_002077035.3_NA12878_prelim_3.0_genomic.fna.gz gunzip GCA_002077035.3_NA12878_prelim_3.0_genomic.fna.gz mv GCA_002077035.3_NA12878_prelim_3.0_genomic.fna NA12878.fa 第二步: 在https://bionanogenomics.com/support/software-downloads/下载Solve软件, 服务器要满足如下需求: Python=2.7.x perl=5.14.x或5.16.x R > 3.1.2,并且安装data.table, igraph, intervals, MASS, parallel, XML, argparser glibc >= 2.14 和 gcc库 至少有一个256G节点的内存,最好有一些32G内存的小节点 tar -zxvf Solve3.3_10252018.tar.gz 解压缩后里有如下几个文件夹 cohortQC: 主要是MQR运行脚本 HybridScaffold: 单酶系统和双酶系统混合组装工具脚本 Pipeline:从头组装的脚本 RefAligner:用于序列联配和组装 RefGenome:hg19和hg38的cmp文件 SVMerge: 用于合并单酶系统得到SV结果 UTIL: 运行从头组装的一些实用shell脚本,可以根据需要进行修改 VariantAnnotation: 对找到的SV进行注释 VCFConverter: 将SMAP和SVMerge的结果输出成VCF格式 目前的主流Bionano设备已经是Sapjyr,BNX文件产生于Saphyr,经由Bionano Access 下载到本地。 从公司拿到的是"RawMolecules.bnx"文件, 不过我们练习用的数据是"output/all.bnx.gz", mkdir test mv output/all.bnx.gz test zcat all.bnx.gz| head -n 20000 | grep "# Run Data" | wc -l # 281 我们发现发现数据集来自于281个通道。 Label SNR 过滤: 过滤信噪比较低的分子,信噪比低意味着质量差。 有如下几个情况,不需要做Label SNR 过滤,或在你BNX文件的"#rh"部分有"SNRFilterType"定义,就不需要过滤 人类样本不需要过滤。 Bionanao Access 1.2 以后新的图像检测算法得到的BNX文件不需要SNR过滤. 对于AutoDetect或Irysview处理过的数据,默认会进行label SNR过滤,处理之后就是Molecules.bnx 由于我们是人类数据集,因此下面的代码就不需要运行了,并且绝大部分情况下也用到下面的命令。 perl /opt/biosoft/Solve3.3_10252018/cohortQC/10252018/filter_SNR_dynamic.pl -i RawMolecules.bnx -o Molecules_filter.bnx -P diag_hist.pdf > snr.log & 分子长度过滤: 过滤短与某个阈值的分子,公司一般会只保留100kb或120Kb以上的分子(取决于数据量,数据越多,阈值越高) gunzip all.bnx.gz RefAlinger -i all.bnx -minlen 120 -merge -o output -bnx > run.log & 在输出的内容中,注意如下部分 Final maps=1147524, sites=46764680, length= 268845841.028kb (avg= 234.283 kb, label density= 17.395 /100kb) 表明过滤后,还有113万条分子,涉及到4676万标记,总测序量为258G,标记密度是17.395/100kb,平均分子长度大于234Kb.. 测序深度等于总测序量除以基因组大小,这是人类基因组,按照3G计算,那么深度就是80X. 过滤后平均分子长度应大于200Kb。 标记密度不能过高,过高会因分辨率不够而无法区分,过低则无法用于比对。一般DLS在10~25 , NRLS在8~15. 组装评估: 在正式组装之前,我们还需要判断下当前数据是否满足组装要求。 为了获取所需的评估参数,得将过滤后的BNX文件和基因组模拟模切得到的CMAP进行比对 第一步: 对基因组序列模拟酶切,得到CMAP文件 perl /opt/biosoft/Solve3.3_10252018/Pipeline/10252018/fa2cmap_multi_color.pl -i ../NA12878.fa -e cttaag 1 mv ../NA12878_CTTAAG_0kb_0labels.cmap . CMAP的格式比较简单,说明如下: 第二步: 用align_bnx_to_cmap.py进行比对。 bionano光学图谱比对的基本原理是基于标记的相对位置。 python /opt/biosoft/Solve3.3_10252018/Pipeline/10252018/align_bnx_to_cmap.py \ --prefix human \ --mol molecules120k.bnx \ --ref NA12878_CTTAAG_0kb_0labels.cmap \ --ra /opt/biosoft/Solve3.3_10252018/RefAligner/7915.7989rel \ --nthreads 80 \ --output prealign \ --snrFilter 2 \ --color 1 参数说明可自行阅读帮助说明。我们重点关注输出结果中如下几方面内容: “contigs/alignmolvref/alignmol_log_simple.txt”里的“Fraction of aligned molecules”和"Effective coverage of reference (X)". 我们要判断数据是否符合最低的比对率。比对率和基因组实际情况有关(组装质量,错误率,重复坍缩情况)。对于人类基因可以达到90%以上,对于不怎么完整度的基因组,即便Bionano的质量很高,比对率可能也只有30%~40%(仅统计150 kb 的分子) "alignments.tar.gz", 里面包含的三个文件可以输入到Bionano提供的另一个工具Access中进行可视化,注意导入要选择"Anchor to Molecules"。 由于人类基因组足够的大,因此需要等待一段时间才能处理完成,之后就可以对比对结果有一个更加直观的了解。 那么合格后的数据应该如何组装呢?请等待后续教程!

常用命令: sort学习笔记

本文的sort命令是GNU版本(8.22), 和BSD的sort不同 sort是我最常用Linux命令之一,它的功能就是排序,一般后面还会和uniq搭配,对数据进行去重。 下面的操作假设你有一个文件,叫做chr.txt, 内容如下, 不同列之间用制表符分隔 Chr3 20251812 20254323 + Chr1 471971 473336 - Chr3 21701520 21703114 + Chr3 18709613 18710836 + Chr5 25645209 25646845 - Chr1 3055231 3056997 - Chr1 28881539 28885023 + 最简单的用法就是只给一个文件作为输入,默认就是从左到右逐个字符的进行比较, $ sort chr.txt Chr1 28881539 28885023 + Chr1 3055231 3056997 - Chr1 471971 473336 - Chr3 18709613 18710836 + Chr3 20251812 20254323 + Chr3 21701520 21703114 + Chr5 25645209 25646845 - 默认是从小到大,也可以用-r, --reverse实现从大到小进行排序 $ sort -r chr.txt Chr5 25645209 25646845 - Chr3 21701520 21703114 + Chr3 20251812 20254323 + Chr3 18709613 18710836 + Chr1 471971 473336 - Chr1 3055231 3056997 - Chr1 28881539 28885023 + 如果你想指定某一列进行排序,需要用到参数-k, --key=KEYDEF, 其中KEYDEF就是你指定的列 $ sort -k2 chr.txt Chr3 18709613 18710836 + Chr3 20251812 20254323 + Chr3 21701520 21703114 + Chr5 25645209 25646845 - Chr1 28881539 28885023 + Chr1 3055231 3056997 - Chr1 471971 473336 - 你会发现一个不对劲的情况,为啥28881539会比3055231小?因为默认情况下,sort将数字当做字符处理,从左到右,逐个比较字符,在第一个数字上,2是比3小,因此28881539在前3055231在后。如果你需要按照数值大小进行排序,那么就需要一个额外参数-n,--numeric-sort $ sort -nk2 chr.txt Chr1 471971 473336 - Chr1 3055231 3056997 - Chr3 18709613 18710836 + Chr3 20251812 20254323 + Chr3 21701520 21703114 + Chr5 25645209 25646845 - Chr1 28881539 28885023 + 那能不能先根据第一列然后根据第二列排序呢? $ sort -k1,1 -k2,2n chr.txt Chr1 471971 473336 - Chr1 3055231 3056997 - Chr1 28881539 28885023 + Chr3 18709613 18710836 + Chr3 20251812 20254323 + Chr3 21701520 21703114 + Chr5 25645209 25646845 - 你可能对这里面 -k1,1 -k2,2n感觉特别的奇怪,这是因为我还没有说明KEYDEF到底是什么 KEYDEF的定义是 F.C,F[OPTS]]。这里[]就是可填可不填的意思,例如一开始的k2就是只填了第一个F。F表示filed, 也就是列,这里一共有两个F,分别用于排序键值的起始位置列和结束位置列。举个例子 $ sort -k1,2 chr.txt Chr1 28881539 28885023 + Chr1 3055231 3056997 - Chr1 471971 473336 - Chr3 18709613 18710836 + Chr3 20251812 20254323 + Chr3 21701520 21703114 + Chr5 25645209 25646845 - 效果就是将第一列和第二列进行合并,然后从左到右进行逐个比较。 .C适用于Chr23 和 Chr2这种字符和数字混合的情况,我们希望Chr2 < Chr23,但是默认是Chr23 > Chr2. .C表示从第C个位置开始,当做数字处理 $ cat test.txt Chr23 $ sort -k1.4 test.txt Chr23 最后OPTS表示该列的排序类型,可选项是 [bdfgiMhnRrV], 我比较常用的就是n和r。 那么 -k1,1 -k2,2n就很好理解了, 就是先按照第一列排序,然后按照第二列排序,其中第二列是数值型数据。 除了按数字排序,sort还支持按照月份缩写-M,--month-sort, 例如'JAN'<'DEC', 人类可读数值-h,--human-numeric-sort, 例如 2K < 1G。 如果你的输入内容是 0.04,1,123,1e-10,3.23 这种数字记录类型,sort也提供-g,--general-numeric-sort,对其排序 $ sort numbers.txt 1e-10 $ sort -g numbers.txt 1e-10 还有一个和软件版本开发相关的排序-V, --version-sort。你是不是觉得1.21 的版本比 1.3新。然而在软件开发领域,1.3 比1.21新 $ sort -V version.txt 其他你可能会用得到的参数 -b: 获取开头的空白 -f: 忽略大小写 --parallel=N: 并行 -o, --output=FILE: 指定输出文件,而非标准输出 -c, --check: 检查是否排序, 输出第一个排序不对的位置 -C, --check=quiet, --check=silent: 检查是否排序,没有任何输出,echo $?返回1 -t, --field-separator=SEP: 默认是空格分隔,可以指定分隔符 -m,--merge: 合并已经排序的文件

如何在shell脚本中控制任务投递

如果只有一个样本,或者样本量不大的情况下,我会选择一次性投递所有的任务。但是如果有100个以上的样本,那我就得谨慎考虑。 用 snakemake 很好解决这个问题,它会按照你给定的任务数和CPU数,确定每次投递多少任务。但是,有些时候任务比较简单,比如说gzip压缩,我就不想写一个snakemake脚本,就想用shell完成。 解决这个问题的关键命令是wait, 它的作用是等待当前的任务完成,和条件语句进行搭配的话,就能够实现每次只运行一定量的任务。举个例子, 我想用gzip对这些文件进行压缩 直接用下面的命令,可能会影响到其他用户 ls *.fastq | while read fastq ; do echo "gzip $fastq &"; done 这里没有实际运行,只是用echo来示意。 而比较合适的做法是下面这种 #!/bin/bash count=0 ls *.fastq | while read fastq gzip $fastq & count=$((count+1)) if [ $((count % 4 )) -eq 0 ]; count=0 每次循环时,会把任务会放入后台,之后检查当前的计数是不是4的倍数。如果是的话,等任务完成,在运行后续的任务,不是的后直接运行下一个任务。 这个小技巧要感谢 李广伟师兄。

HaploMerger2: 从高杂合二倍体基因组组装中重建单倍型

本文只是按照自己的需求翻译了HaploMerger2提供的手册部分内容。HaploMerger2的帮助文档写的非常好,一定要花点时间去读啊! HaploMerger2的分析流程如下 重建单倍体组装中的等位基因关系 检测并纠正二倍体组装中的错连(mis-join) 重建2个单倍型组装 进一步对单倍型组装进行桥接 检测并移除串联错配 图片流程: 对于如何使用HaploMerger2获取高质量分型单倍型组装,作者从Canu和Falcon/Falcon-unzip的结果分别给了建议 如果是Canu,尽量避免基因组坍缩,尽可能获取完整的分离的二倍型组装;然后用HaploMerger2获取高质量的参考基因组和可选单倍型组装;然后用HapCUT2或者其他分型工具基于参考单倍型组装去获取高质量单倍型组装, 对于Falcon结果,将Falcon的输出结果p-tigs和a-tigs进行合并,然后将结果传给HaploMerger2 在https://github.com/mapleforest/HaploMerger2/releases下载HaploMerger2,解压缩后会得到很多的文件夹 bin: HaploMerger2的核心程序 chainNet gapCloser lastz SSPACE-STANDARD winMasker project_template: 配置文件 project_template里的运行脚本用的都是相对路径,也就是../bin/软件名, 这其实很不方便进行脚本调用,因此,推荐用下面的方法将../bin替换成你的实际的bin路径。(如果担心操作失误,可以备份project_template) cp -r project_template project_template_bck cd project_template # 确认要被更改的信息 grep '\.\.' hm* # 用sed进行更改 sed -i 's/\.\./\/opt\/biosoft\/HaploMerger2/' hm* 注: 我的路径是/opt/biosoft/HaploMerger2,你需要按照自己的情况进行调整。 最后还需要保证lastz, chainNet, winMasker, SSPACE, GapCloser这些软件应该在环境变量中。 最低要求: 基因组压缩后的文件。 因为我没有处理过10k, 20k, 50k文库用于scaffold的项目,而且D流程会删除基因组序列, 所以本次实战只运行A, B这2个流程. 在HaploMerger2目录下新建一个项目文件,命名为example,之后下载分析所用的数据。然后将你的基因组序列移动到example里,例如 contig.fa mkdir project_your_name cd project_your_name cp /path/to/your/contig.fa . 从project_template拷贝流程对应的脚本和配置文件 cp ../project_template/hm.batchA* . cp ../project_template/hm.batchB* . cp ../project_template/hm.batchD* . cp ../project_template/all_lastz.ctl . cp ../project_template/scoreMatrix.q . 第一步:对重复序列进行软屏蔽 这是必须要做的一步,能提高后续分析的准确性。推荐用windowmasker。 构建能重复使用的屏蔽文库 windowmasker \ -checkdup true \ -mk_counts \ -in contigs.fa \ -out masking-library.ustat \ -mem 6500 然后用构建的文库对基因组进行重复序列软屏蔽,也就是将重复序列从大写改成小写。 windowmasker \ -ustat masking-library.ustat \ -in contigs.fa \ -out contigs_wm.fa \ -outfmt fasta \ -dust true 然后将输出文件进行压缩 gzip contigs_wm.fa 第二步: 移除组装结果中主要的错连 如果你发现运行脚本的时候只用了几秒钟就结束了,那么基本上就是你的环境没有配置好。 hm.batchA1是对文件进行拆分,然后用LASTZ进行all-vs-all的全基因组比对。比对结果会存放在contigs_wm.contigs_wmx.result/raw.axt(contigs_wm是前缀), 根据你设置的特异性不同,占用的空间大小也不同,在运行完第二个脚本(A2),可以删除。 sh ./hm.batchA1.initiation_and_all_lastz contigs_wm hm.batchA2是利用ChainNet算法创建互为最优的单覆盖(reciprocally-best single-coverage)全基因组联配 sh ./hm.batchA2.chainNet_and_netToMaf contigs_wm 这两步都可以通过修改其中的threads参数来提高运行速度,其中LASTZ每个线程要求约1G内存,而chainNet需要至少4-8G内存。 如果为了降低非同源联配,可以修改hm.batchA1中的identity参数,默认是80.如果你的基因组杂合度也就是4~5%左右,那么设置88%~92%能够获得更加严格的过滤。 这两步有两个非常重要的配置文件,all_lastz.ctl和scoreMatrix.q。其中all_lastz.ctl控制LASTZ,例如你可以将--step=1改成--step=20来提高LASTZ的运行速度。默认的参数其实已经挺不错了,如果想要修改的话,建议先去看看LASTZ的手册。 scoreMatrix.q里面记录的是核酸替换的得分矩阵,用于LASTZ和chainNet中。默认的得分矩阵来自于Chinese amphioxus的二倍体组装(GC 41%, 杂合度4~5%)。为了更好的区分等位基因和非等位基因的联配,建议参考后面的如何自定义LASTZ得分矩阵部分进行构建。 hm.batchA3会调用5个Perl程序,完成5个任务 预过滤互为最优的全基因组chainNet联配 将全基因组联配转成有向图(DGA graph 使用贪婪算法遍历有向图,寻找错连 将序列在错连的位置进行打断 基于有向图输出二倍体组装 sh hm.batchA3.misjoin_processing contigs_wm 这一步可以进行调整的参数如下 scoreSchme: 计分规则,目前使用默认的score即可,分数越高说明联配质量越高 filterScore_1, escapeFilter, filterScore_2和NsLCsFilter: 大部分低得分的联配可能都是假的,应该被归为噪音。有三个参数用于进行控制过滤规则。 对于第二步,作者建议迭代2~3次,保证尽可能的找到错连的部分。同时还建议尝试下调aliFilter和overhangFilter(从50kb到40kb即可,低于30kb要十分小心) 假如你迭代了三次,那么每次的结果应该是contigs_wm > contigs_wm_A > contigs_wm_A_A > contigs_wm_A_A_A 第三步: 创建二倍型组装 先运行hm.batchB1.initiation_and_all_lastz和hm.batchB2.chainNet_and_netToMaf sh hm.batchB1.initiation_and_all_lastz contigs_wm_A_A_A sh hm.batchB2.chainNet_and_netToMaf contigs_wm_A_A_A 前两个脚本和之前的batchA的前两个脚本类似,除了两点: 二倍体输入文件应该是第一步的输出结果 chaiNet联配会输出MAF文件,可以用Gmaj查看 运行 hm.batchB3.haplomerger sh hm.batchB3.haplomerger contigs_wm_A_A_A 第三个脚本做了以下四件事情 预过滤互为最优的全基因组chainNet联配 将全基因组联配转成有向图(DGA graph) 使用贪婪算法遍历和线性化DGA 基于线性化的DGA输出参考单倍型和可选单倍型 hm.batchB3和hm.batchA3的参数类似。作者发现,提高filterScore_2 和escapeFilter能够避免将原本不是等位基因的两条序列进行合并,但是却会导致单倍型组装结果丢失正确的联配信息,从而导致包含了更多的冗余序列。部分的联配中主要是N和重复序列,使用NsLCsFilter能够过滤这些联配。 同时HaploMerger还能够将2个存在重叠区域(overlapping region)的序列进行合并, 参数minOverlap控制合并操作的最小联配长度(不含N, gaps和InDels)。默认设置是0,也就是将任何可能合并的序列都进行合并。由于低分联配已经被filterScore_2过滤,所以最小联配长度实际上是大于filterScore_2。提高minOverlap会让结果更加可信,不过肯定也会损失可能的连接。 因为任何合并都会将两个单倍型合并,或者你不想混合单倍型或在单倍型间产生交换(generate switches between haplotypes),那么你可以设置minOverlap= 99999999 这一步会在contigs_wm_A_A_A.contigs_wm_A_A_Ax.result文件下输出三个重要文件,分别是 optiNewScaffolds.fa.gz: 有联配支持的新参考序列 optiNewScaffolds_alt.fa.gz: 有联配支持的可选参考序列 unpaired.fa.gz: 没有联配支持的序列 第四个脚本是优化没有联配支持的序列 sh hm.batchB4.refine_unpaired_sequences contigs_wm_A_A_A unpaired.fa.gz里存放的是没有联配支持的序列,有以下几个原因会导致该情况 一些序列在二倍型组装中没有对应的等位序列 等位序列/单倍型坍缩成单个序列 一些真实的联配由于信号太弱被过滤 一些联配包含N和重复序列,导致它被过滤 一些联配和串联重复,倒置和易位,所以在线性化的DGA中被忽略 在图遍历是,部分序列会被过滤(例如序列中的无法联配自由末端) 考虑到这些情况,我们会认为这些序列中大部分是冗余序列。所以bactchB4脚本就是找到这些冗余序列,然后进行过滤。该脚本会将unpaired.fa.gz和optiNewScaffolds.fa.gz进行比对,然后移除其中潜在冗余序列。这会得到新的fasta文件,即unpaired_refined.fa.gz。几个可供调整的参数 threads: 控制线程数 identity: 控制特异性和敏感性 maskFilter: 过滤只有N和重复序列的scaffold redudantFilter: 过滤在optiNewScaffolds.fa.gz有对应等位序列的scaffold 最后一步就是合并。 将optiNewScaffolds.fa.gz 与 unpaired_refined.fa.gz合并得到参考单倍型contigs_wm_A_A_A_ref.fa.gz, 将optiNewScaffolds_alt.fa.gz 与 unpaired_refined.fa.gz合并得到contigs_wm_A_A_A_alt.fa.gz sh hm.batchB5.merge_paired_and_unpaired_sequences contigs_wm_A_A_A 最后的结果是contigs_wm_A_A_A_ref.fa.gz和contigs_wm_A_A_A_alt.fa.gz 如何自定义LASTZ得分矩阵 作者为了方便我们推断LASTZ_D的得分矩阵,封装了一个脚本 lastz_D_Wrapper.pl 我们需要根据基因组大小将基因组序列分成2个部分,一部分只有5~15%的序列,而另一部分则是剩下的序列。一般而言,你应该选择最长的序列用于第一部分。 lastz_D_Wrapper.pl --target=part1.fa.gz --query=part2.fa.gz --identity=90 这里的参数主要是--identity=90, 表示lastz_D只会使用大于90%相似度的联配,该值越高表示严格度越高。作者推荐当杂合率高在4%~5%的时候选择90%,在1%~3%的时候设置为94%~96%。 运行完之后,将 part1.part2.raw.time.xx_xx.q的内容复制到scoreMatrix.q中。 注1: lastz_D速度很慢而且不支持并行,因此只要用10%以内的序列作为target就行。序列多了反而找到等位基因的概率低了。 注2: 有些时候lastz_D还会诡异的停住,然后输出奇怪的结果。此外,不同的序列还会有不同的推断结果。因此,part1.fa.gz可能要选择不同的序列,然后运行程序。之后,你可以从不同的得分矩阵中筛选结果。 注3: 这一步和基因组复杂度密切相关,杂合度越高,运行时间越长。 如果在分析中用到了HaploMerger2,请引用下面的文章 Shengfeng Huang, et al. HaploMerger2: rebuilding both haploid sub-assemblies from high-heterozygosity diploid genome assembly. Bioinformatics. 2017. Shengfeng Huang, et al. HaploMerger: reconstructing allelic relationships for polymorphic diploid genome assemblies. Genome Res. 2012, 22(8):1581-1588.

ChIPseeker的upsetplot是怎么写的

在距今2年前的一个日子里,我写了一篇关于upsetR的简单教程。两年过去了,这篇教程是百度搜索的第一位。 在这篇教程中,我定了一个计划 下一期写一篇Y叔的upsetplot是如何写的。 然而两年过去了,我说的那个教程都没有出现,今天有人提问upsetR的时候,我又看到了这个计划,羞愧难当,遂有了这篇教程。 通常而言,看一个函数的源代码的方式就是直接在R中输入这个函数名 > upsetplot standardGeneric for "upsetplot" defined from package "enrichplot" function (x, ...) standardGeneric("upsetplot") <bytecode: 0x0000000047c188c8> <environment: 0x0000000047c10318> Methods may be defined for arguments: x Use showMethods("upsetplot") for currently available ones. 但是对于upsetplot这个函数而言,并不奏效,因为它是一个泛型函数。所谓的泛型函数,指的是对于同一个函数而言,它会根据不同数据对象调取实际的函数。例如print函数,既可以用来输出字符串,还可以输出图片。 那么如何查看它的实际函数呢?我们可以根据提示,showMethods("upsetplot")去找它到底有哪些可用函数 > showMethods("upsetplot") Function: upsetplot (package enrichplot) x="csAnno" x="enrichResult" 不难发现,它提供了两种数据类型的展示方法,一种是csAnno,另一种是enrichResult.对于csAnno,我们可以用包名:::不可见函数的形式查看到实际代码 ChIPseeker:::upsetplot.csAnno function (x, sets = NULL, order.by = "freq", sets.bar.color = NULL, vennpie = FALSE, vp = viewport(x = 0.6, y = 0.7, width = 0.6, height = 0.6), ...) y <- x@detailGenomicAnnotation y <- as.matrix(y) y[y] <- 1 y <- as.data.frame(y) if (is.null(sets)) { sets <- c("distal_intergenic", "downstream", "threeUTR", "fiveUTR", "Intron", "Exon", "Promoter") if (vennpie && is.null(sets.bar.color)) { sets.bar.color <- c("#d95f0e", "#fee0d2", "#98D277", "#6F9E4C", "#fc9272", "#9ecae1", "#ffeda0") if (is.null(sets.bar.color)) { sets.bar.color <- "black" if (vennpie) { plot.new() upset(y, sets = sets, sets.bar.color = sets.bar.color, order.by = order.by, ...) pushViewport(vp) par(plt = gridPLT(), new = TRUE) vennpie(x) popViewport() else { upset(y, sets = sets, sets.bar.color = sets.bar.color, order.by = order.by, ...) 为了能够解释代码,我们需要先运行测试数据,获取能够进行演示的数据 require(TxDb.Hsapiens.UCSC.hg19.knownGene) txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene peakfile <- system.file("extdata", "sample_peaks.txt", package="ChIPseeker") peakAnno <- annotatePeak(peakfile, tssRegion=c(-3000, 3000), TxDb=txdb) x <- peakAnno # x是csAnno类 y <- x@detailGenomicAnnotation用来提取详细的注释信息 y <- as.matrix(y)转换成矩阵,方便后面的y[y]<-1将所有为TRUE的结果转成1,在y<- as.data.frame(y)的时候让FALSE转成了0。 如果是我的话,我应该会用apply y <- x@detailGenomicAnnotation y2 <- apply(y, 2, as.integer) y2 <- as.data.frame(y2) 这里的y的数据结构刚好是的upset需要的输入,所以假如Y叔没有提供upsetplot函数的话,那么代码你可以直接用upset(y)作图 upset(y) 下面代码定义需要提供给upset函数进行集合运算的属性,默认就是ChIPseeker注释的8个属性 之后的代码是根据是否绘制饼图走向两条分支。 在绘制饼图的循环里,Y叔先用graphics的plot.new新建了一个绘图框架,方便在这个绘图框架中作图。 plot.new() 再用upset把图画好 # 绘图的准备工作 sets <- c("distal_intergenic", "downstream", "threeUTR", "fiveUTR", "Intron", "Exon", "Promoter") sets.bar.color <- c("#d95f0e", "#fee0d2", "#98D277", "#6F9E4C", "#fc9272", "#9ecae1", "#ffeda0") sets.bar.color <- "black" order.by <- "freq" library(UpSetR) upset(y, sets = sets, sets.bar.color = sets.bar.color, order.by = order.by) 之后的pushViewport(vp)会调出一个视图对象,也就是vp = viewport(x = 0.6, y = 0.7, width = 0.6, height = 0.6), 指定了后续要画图的位置和大小. par(plt = gridPLT(), new = TRUE)设置了绘图的基本参数。gridPLT函数来自于gridBase,他的作用就是返回用于设置plt的值。 最终用vennpie(x)绘制饼图,然后用popViewport()回到原来视图,结束画图。

eggnog-mapper实现功能注释 eggNOG-Mapper介绍 通常功能注释的思路都是基于序列相似性找直系同源基因,常见的方法就是BLAST+BLAST2GO, 或者是InterProScan。eggNOG-mapper的作者认为这种方法不够可靠,毕竟你有可能找到的的是旁系同源基因。近期对多个工具的整体评估发现eggNOG(evolutionary genealogy of genes: Non-supervised Orthologous Groups)在区分旁系同源基因和直系同源基因上表现不错,因此基于eggNOG数据库开发了eggNOG-mapper工具,用于对新序列进行功能注释。 eggNOG-mapper的算法实现如下: 第一步:序列比对。首先,每条蛋白序列用HMMER3在整理的eggNOG数据库中搜索。由于每个HMM匹配都和一个功能注释的eggNOG OG对应,这一步就提供了初步的注释信息。之后,每条蛋白序列用phmmer在最佳匹配的HMM对应的一组eggNOG蛋白中进一步搜索。最后,每条序列的最佳匹配结果以 seed ortholog 形式存放,用于获取其他直系同源基因。目前eggNOG HMM数据库中拥有1,911,745个OG,覆盖了1,678种细菌,115种古细菌,238种真核物种以及352种病毒。除了HMMER3外,还而可用DIAMOND直接对所有的eggNOG蛋白序列进行搜索,它的速度更快,适合类似于宏基因组这类大数据集,或者是已有物种和eggNOG所收集的物种比较近。当然服务器性能强大的话,还是有限选择HMMER3. 第二步:推测直系同源基因。每个用于检索的蛋白序列的最佳匹配序列会对应eggNOG的一个蛋白, 这些蛋白基于预分析的eggNOG进化树数据库会提取一组更加精细的直系同源基因。这一步还会根据bit-screo或E-value对结果进行一次过来,剔除同源性不高的结果 第三步:功能注释。用于搜索的蛋白序列对应的直系同源基因的功能描述就是最终的注释结果。比如说GO, KEGG, COG等。 安装本体eggnog-mapper之前,需要先保证服务器上已经安装Python3.7(BioPython模块), wget, HMMER3 和/或 DIAMON,此外你还得保证70G的空间用于存放注释数据库和FASTA文件,对于真核生物至少保证90G的服务器内存。之后从GitHub上将软件下载到本地 git clone https://github.com/jhcepas/eggnog-mapper.git 之后需要下载所需要的数据库。eggNOG提供了107个不同物种的HMM数据库(xxxNOG)以及三个优化数据库, euk对应真核生物,bact对应细菌, arch对应古细菌(Archeabacteria), 以及一个病毒数据库(viruses). 这三个优化数据库包含了属于该分类内的所有物种的HMM模型。 cd eggnog-mapper ./download_eggnog_data.py euk 此处下载的是真核生物 。下载过程中,它会反复询问你是否要下载某一类数据,我一律选择是。如果你的服务器有多个盘,安装软件的分区不够大,可以在数据盘中创建文件夹进行软链接 cd eggnog-mapper mkdir /data/database/EGGNOG-DB rmdir data && ln -s /data/database/EGGNOG-DB data eggnog-mapper用起来非常的简单,你需要提供蛋白序列作为输入 #假如我们现在都仍在软件安装的路径下 python emapper.py -i test/p53.fa --output p53_maNOG -d euk python emapper.py -i test/p53.fa --output p53_maNOG -d maNOG python emapper.py -i test/p53.fa --output p53_maNOG -d maNOG --usemem --cpu 10 eggnog-mapper默认是以HMMER进行序列搜索,尽管这可以通过-m diamond更改成DIAMOND,但结果中就会缺少一些信息列。HMM的搜索参数,有--hmm_maxhits, --hmm_evalue, --hmm_score, --hmm_qcov和--Z, 一般用默认就行了。 -d/--database表示搜索数据库的类型,既可以是"euk, bact, arch" 三大类的其中一种,也可以是"eggnog-mapper/data/OG_fasta"中的小类,例如"maNOG"表示的就是哺乳动物的NOG。此外也可以是自定义HMM数据库, -d /path/to/pfam.hmm 如果服务器内存比较大,线程比较多,可以用--usemem和 --cpu 线程数提高运行速度。 --output表示输出文件的前缀,默认输出在当前文件夹下,--output_dir可以更改为其他文件路径。--resume表示任务重启后可以跳过之前已经完成的部分, 而--override则表示覆盖原先的输出结果。 如果你计划使用同一个数据库同时对多个物种进行注释(大部分人都没有这种需求),那么你需要以服务器模式启动eggnog-mapper,然后就可以同时启动多个任务, # 终端1 python emapper.py -d arch --cpu 10 --servermode # 新建一个终端 python emapper.py -d arch:localhost:51600 -i test/polb.fa -o polb_arch eggnog-mapper会生成三个文件, [project_name].emapper.hmm_hits: 记录每个用于搜索序列对应的所有的显著性的eggNOG Orthologous Groups(OG). 所有标记为"-"则表明该序列未找到可能的OG [project_name].emapper.seed_orthologs: 记录每个用于搜索序列对的的最佳的OG,也就是[project_name].emapper.hmm_hits里选择得分最高的结果。之后会从eggNOG中提取更精细的直系同源关系(orthology relationships) [project_name].emapper.annotations: 该文件提供了最终的注释结果。大部分需要的内容都可以通过写脚本从从提取,一共有13列。 [project_name].emapper.annotations每一列对应的记录如下: query_name: 检索的基因名或者其他ID sedd_eggNOG_ortholog: eggNOG中最佳的蛋白匹配 seed_orholog_evalue: 最佳匹配的e-value seed_ortolog_evalu: 最佳匹配的bit-score predicted_gene_name: 预测的基因名,特别指的是类似AP2有一定含义的基因名,而不是AT2G17950这类编号 GO_term: 推测的GO的词条, 未必最新 KEGG_KO: 推测的KEGG KO词条, 未必最新 BiGG_Reactions: BiGG代谢反应的预测结果 Annotation_tax_scope: 对该序列在分类范围的注释 Matching_OGs: 匹配的eggNOG Orthologous Groups best_OG|evalue|score: 最佳匹配的OG(HMM模式才有) COG functional categories: 从最佳匹配的OG中推测出的COG功能分类 eggNOG_HMM_model_annotation: 从最佳匹配的OG中推测出eggNOG功能描述 如果打算做富集分析,用命令行的cut/awk提取对应的列,过滤掉其中未注释的部分就行了。 eggNOG-Mapper的文章: 如果你专门是做数据分析,用别人的软件的时候最好把文章也看了 官方教程: 简洁明了

使用ALLMAPS进行辅助组装 在从头组装过程中,确定基因组的scaffolds/contig的顺序和朝向是重建染色体非常关键的一步。这一步可以由多种辅助组装策略完成:例如遗传图谱, Hi-C, BioNano光学图谱,10X Chicago 。 一个物种可能会有多个遗传图谱,可以是不同项目中的不同定位群体结果,可以是不同软件如R/QTL, MSTMAP和JOINMAP的分析结果。遗传图谱会因重组率,偏分离(segregation distortion) , PAV(presence-absence variation)和染色体比对多态位点不同而发生变化。每一种图谱都能够提供不同的证据,举个例子,一个scaffold可能在一个图谱中无法被锚定,但是在另一个图谱中可以进行锚定,将这些图谱进行整合就能提最后染色体组装的精确度。 如果只用一个图谱,对scaffold进行排序只是计算量大一点而已,你需要根据图谱中分子标记在每个scaffold的平均距离进行排序就行。ALLMAPS, 正如名字说的那样,就是能够使用所有的图谱证据的工具,它能够计算scaffold的朝向,使其和已有的图谱的共线性关系最大化。他有如下亮点: 可重复性:清晰的可计算目标使得让多种输入图谱的共线性关系能够最大化 灵活性:允许为输入的图谱设置权重,更好的处理冲突 强大性:能使用多种遗传图谱,只需要做最小的转换 易用性:让基因组构建(FASTA和AGP输出)和基因组升级(CHAIN输出)流程化 通常,解决scaffold顺序和朝向(OO)是一种NP问题。因为基因组组装和遗传图谱中都可能存在一些错误,我们的目标只能是找到一个近似解。ALLMAPS将该问题转换成旅行商问题,然后用遗传算法优化scaffold OO. 使用遗传算法优化是为了避免在局部最优上出现瓶颈。遗传图谱最常见的错误是倒置和异位(inversion and translocation) 安装可以用CONDA。 conda install jcvi 先现在一个测试数据集,用于练习了解软件的总体流程。比较尴尬的是,这个数据集放在了Dropbox上,我们都知道这是一个国内不存在的网站,所以我把他转存到了我的腾讯微云上,链接:https://share.weiyun.com/5nwjljN . 第一步:准备输入文件。 首先你得要提供物理图谱和遗传图谱的对应关系,格式为 Scaffold ID, scaffold position, LG, genetic position 你可以在excel表格中进行操作,然后另存为csv为文件,如下所示。你可以发现有一些scaffold里就只有一个遗传标记进行对应,你可以思考下该scaffold最后会如何处理 当然,我们解压缩后的zip文件里就有这些内容,你可以用一些命令工具(less, head, tail)看下数据是如何存放的。 第二步: 将两个图谱进行合并,最后会得到一个权重文件(weights.txt)和输入的bed文件 python -m jcvi.assembly.allmaps merge JMMale.csv JMFemale.csv -o JM-2-test.bed 第三步:对权重文件"weights.txt" 进行调整。weights.txt默认每个输入的图谱的权重都是1。 作者有一个建议就是,你通过检查最后的报告和诊断图,有监督地来重新对每个遗传图谱进行权重赋值。 $ cat weights.txt JMFemale 1 JMMale 1 第四步: 对scaffold进行排序,搭建成准染色体水平 python -m jcvi.assembly.allmaps path -w weights.txt JM-2-test.bed scaffolds.fasta 上面第四步会在当前文件夹下生成许多结果文件, JM-2-test.agp JM-2-test.bed JM-2-test.chain JM-2-test.chr.agp JM-2-test.chr.fasta JM-2-test.fasta JM-2-test.fasta.sizes JM-2-test.lifted.bed JM-2-test.summary.txt JM-2-test.tour JM-2-test.unplaced.agp JM-2-test.unplaced.fasta 可以分为如下几类 首先是组装后的基因组序列 "JM-2-test.fasta" "JM-2-test.agp": 每个scaffold的顺序和朝向,用于上传Genbank "JM-2-test.chain": 用于新旧坐标间的转换,比如说你在之前坐标注释得到GFF文件就可以用lifOver转换到新坐标系下 然后可视化报告。每个染色体都会得到一个对应的pdf文件,可视化展示如下:左图主要关注交叉的线,表示某些marker存在矛盾。右图关注是否有单上升/下降趋势,图中的斜率反应的是物理距离(x轴)相对遗传距离(y轴)的变化,可以认为是重组率。低重组率不容易确定在同一个重组区间内的scaffold的位置和朝向 白色部分为可信区间,灰色部分是存疑区间。 除了默认的输出外,你还可以通过movie得到软件运行过程每次迭代后组装情况,不过要额外安装ffmpeg和parallel python -m jcvi.assembly.allmaps movie -w weights.txt JM-2-test.bed scaffolds.fasta chr23 最后是总结性报告。尽管可视化会给你比较直观的结果,但是具体的参数还是要看JM-2-test.summary.txt. 例如遗传标记的密度,以及有多少序列被锚定到基因组上。 当然ALLMAPS还有一些比较高级的操作: 拆分嵌合contig:嵌合contig指的是一段区域能够比对到多个连锁群中或者染色体中,一个常见的来源就是不同染色体中的重复区域由于过于相似在组装的时候坍缩成了一个。ALLMAPS也提供split进行拆分。 估计gap长度: 默认会用100个N在填充两个scaffold连接的区域,方便Genbank识别其为未知区域。你可以根据遗传距离在不同物理位置上的对应关系,预测不同gap的近似大小然后进行填充,子命令是estimategaps。 在ALLMAPS中使用多种遗传图谱. 光学图谱在着丝粒区间表型不好,主要是里面的串联重复过多,使得限制性内切酶数目减少,遗传图谱中在这个区域依旧有比较多的标记 标记数目和密度。过多计算量大,太小准确性低 基因组组装和建图软件的潜在错误。基因组组装的原始质量对ALLMAPS的影响非常大,尤其是高度片段化的结果可能会导致错误结果。 基因组结构特征也会有影响,例如染色体倒置、异位和片段重复。 参考资料: ALLMAPS: robust scaffold ordering based on multiple maps https://github.com/tanghaibao/jcvi/wiki/ALLMAPS

三代转录组系列:使用Cogent重建基因组编码区

尽管目前已测序的物种已经很多了,但是对于一些动辄几个G的复杂基因组,还没有某个课题组有那么大的经费去测序,所以仍旧缺少高质量的完整基因组,那么这个时候一个高质量的转录组还是能够凑合用的。 二代测序的组装结果只能是差强人意,最好的结果就是不要组装,直把原模原样的转录组给你是最好的。PacBio Iso-Seq 做的就是这件事情。只不过Iso-Seq测得是转录本,由于有些基因存在可变剪切现象,所以所有将一个基因的所有转录本放在一起看,搞出一个尽量完整的编码区。 Cogent能使用高质量的全长转录组测序数据对基因组编码区进行重建,示意图如下: 虽然说软件可以直接通过conda进行安装,但是根据官方文档的流程,感觉还是很麻烦 下面操作默认你安装好了miniconda3, 且miniconda3的位置在家目录下 conda create -n cogent cogent -y source activate cogent conda install biopython pulp conda install networkx==1.10 conda install -c http://conda.anaconda.org/cgat bx-python==0.7.3 pip install -i https://pypi.tuna.tsinghua.edu.cn/simple scikit-image cd ~/miniconda3/envs/cogent git clone https://github.com/Magdoll/Cogent.git cd Cogent git checkout tags/v3.3 git submodule update --init --recursive cd Complete-Striped-Smith-Waterman-Library/src export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$HOME/miniconda3/envs/cogent/Cogent/Complete-Striped-Smith-Waterman-Library/src export PYTHONPATH=$PYTHONPATH:$HOME/miniconda3/envs/cogent/Cogent/Complete-Striped-Smith-Waterman-Library/src cd ../../ python setup.py build python setup.py install 经过上面一波操作或,请用下面这行命令验证是否安装成功 run_mash.py --version run_mash.py Cogent 3.3 此外还需要安装另外两个软件,分别是Minimap2和Mash conda install minimap2 wget https://github.com/marbl/Mash/releases/download/v1.0.2/mash-Linux64-v1.0.2.tar.gz tar xf mash-Linux64-v1.0.2.tar.gz mv mash $HOME/miniconda3/envs/cogent/bin 对待大数据集,你还需要安装Cupcake git clone https://github.com/Magdoll/cDNA_Cupcake.git cd cDNA_Cupcake git checkout -b tofu2 tofu2_v21 python setup.py build python setup.py install 创建伪基因组 让我们先下载测试所用的数据集, mkdir test cd test https://raw.githubusercontent.com/Magdoll/Cogent/master/test_data/test_human.fa ### 第一步: 从数据集中搜索同一个基因簇(gene family)的序列 这一步分为超过20,000 条序列的大数据集和低于20,000条序列这两种情况, 虽然无论那种情况,在这里我们都只用刚下载的测试数据集 第一步:从输入中计算k-mer谱和配对距离 run_mash.py -k 30 --cpus=12 test_human.fa # -k, k-mer大小 # --cpus, 进程数 你一定要保证你的输入是fasta格式,因为该工具目前无法自动判断输入是否是fasta格式,所以当你提供诡异的输入时,它会报错,然后继续在后台折腾。 上面这行命令做的事情是: 将输入的fasta文件进行拆分,你可以用--chunk_size指定每个分块的大小 对于每个分块计算k-mer谱 对于每个配对的分块,计算k-mer相似度(距离) 将分块合并成单个输出文件 第二步:处理距离文件,创建基因簇 process_kmer_to_graph.py test_human.fa test_human.fa.s1000k30.dist test_human human 这里的test_human是输出文件夹,human表示输出文件名前缀。此外如果你有输入文件中每个转录亚型的丰度,那么你可以用-c参数指定该文件。 这一步会得到输出日志test_human.partition.txt,以及test_human下有每个基因family的次文件夹,文件里包含着每个基因簇的相似序列。对于不属于任何基因family的序列,会在日志文件种专门说明,这里是human.partition.txt 如果是超过20,000点大数据集,分析就需要用到Minimap2和Cpucake了。分为如下几个步骤: 使用minimap2对数据进行粗分组 对于每个分组,使用上面提到的精细的寻找family工具 最后将每个分组的结果进行合并 第一步:使用minimap进行分组 run_preCluster.py要求输入文件名为isoseq_flnc.fasta, 所以需要先进行软连接 ln -s test_human.fa isoseq_flnc.fasta run_preCluster.py --cpus=20 为每个分组构建基因簇寻找命令 generate_batch_cmd_for_Cogent_family_finding.py --cpus=12 --cmd_filename=cmd preCluster.cluster_info.csv preCluster_out test_human 得到的是 cmd 文件,这个cmd可以直接用bash cmd运行,也可以投递到任务调度系统。 最后将结果进行合并 printf "Partition\tSize\tMembers\n" > final.partition.txt ls preCluster_out/*/*partition.txt | xargs -n1 -i sed '1d; $d' {} | cat >> final.partition.txt 第二步:重建编码基因组 上一步得到每个基因簇, 可以分别重构编码基因组,所用的命令是reconstruct_contig.py 使用方法: reconstruct_contig.py [-h] [-k KMER_SIZE] [-p OUTPUT_PREFIX] [-G GENOME_FASTA_MMI] [-S SPECIES_NAME] dirname 如果你手头有一个质量不是很高的基因组,可以使用-G GENOME_FASTA_MMI和-S SPECIES_NAME提供参考基因组的信息。毕竟有一些内含子是所有转录本都缺少的,提供基因组信息,可以补充这部分缺失。 由于有多个基因簇,所以还需要批量运行命令 reconstruct_contig.py generate_batch_cmd_for_Cogent_reconstruction.py test_human > batch_recont.sh bash batch_recont.sh 第三步: 创建伪基因组 首先获取未分配的序列, 这里用到get_seqs_from_list.py脚本来自于Cupcake, 你需要将cDNA_Cupcake/sequence添加到的环境变量PATH中。 tail -n 1 human.partition.txt | tr ',' '\n' > unassigned.list get_seqs_from_list.py hq_isoforms.fasta unassigned.list > unassigned.fasta 这里的测试数据集里并没有未分配的序列,所以上面这一步可省去。 然后将未分配的序列和Cogent contigs进行合并 mkdir collected cd collected cat ../test_human/*/cogent2.renamed.fasta ../unassigned.fasta > cogent.fake_genome.fasta 最后,将冗余的转录亚型进行合并。 这一步我们需要用Minimap2或者GMAP将Iso-seq分许得到的hq_isoforms.fasta回帖到我们的伪基因组上。 关于参数选择,请阅读Best practice for aligning Iso Seq to reference genome: minimap2, GMAP, STAR, BLAT ln -s ../test_human.fa hq_isoforms.fasta minimap2 -ax splice -t 30 -uf --secondary=no \ cogent.fake_genome.fasta hq_isoforms.fasta > \ hq_isoforms.fasta.sam 然后可以根据collapse tutorial from Cupcake将冗余的转录亚型合并 sort -k 3,3 -k 4,4n hq_isoforms.fasta.sam > hq_isoforms.fasta.sorted.sam collapse_isoforms_by_sam.py --input hq_isoforms.fasta -s hq_isoforms.fasta.sorted.sam \ -c 0.95 -i 0.85 --dun-merge-5-shorter \ -o hq_isoforms.fasta.no5merge 由于这里没有每个转录亚型的丰度文件cluster_report.csv,所以下面的命令不用运行, 最终结果就是hq_isoforms.fasta.no5merge.collapsed.rep.fa get_abundance_post_collapse.py hq_isoforms.fasta.no5merge.collapsed cluster_report.csv filter_away_subset.py hq_isoforms.fasta.no5merge.collapsed 如果运行了上面这行代码,那么最终文件就应是hq_isoforms.fasta.no5merge.collapsed.filtered.*

使用Trinity进行转录组组装

Trinity是Broad Institute和Hebrew University of Jerusalem开发的RNA-Seq数据 转录组组装工具,包括三个模块, Inchworn(尺蠖): 将RNA-seq数据组装成单个转录本,通常是主要转录亚型的全长转录本 Chrysalis(蛹): 这一步将上一步得到contig进行聚类,对于每个聚类构建完整的德布罗意图(_de Bruijin_ graph)。每个转录本表示的是给定基因或者一组有着共同序列的基因的全部转录组成。 之后会根据图中不相交的点对全部短读数据进行拆分 Butterfly(蝴蝶): 并行处理各个图(graph), 追踪每个图中的短读和配对短读的路径,最后报告可变剪切亚型的全长转录本,并且区分出旁系同源基因的转录本 如果不能理解上面这段话,就尝试理解下面这张图吧 当然如果示意图也让你不好理解的话,那就直接用软件吧,反正这些流程图的目标就是想告诉你,“用我,没毛病” 软件安装用bioconda就行了。 conda create -n Trinity trinity -y source activate Trinity 当你在命令行敲出Trinity后,他就会输出一大堆信息。那么多信息分成3个部分: 必须参数:包括--seqType表示输入序列类型,--max_memory允许使用最大内存(一般64G),还有输入数据的所在位置 可选参数:包括链特异性测序参数--SS_lib_type, 线程数--CPU, 允许的最低组装contig长度--min_contig_length, 是否标准化--no_normalize_reads等 常见用法说明 Trinity --seqType fq --max_memory 50G \ --left condA_1.fq.gz,condB_1.fq.gz,condC_1.fq.gz \ --right condA_2.fq.gz,condB_2.fq.gz,condC_2.fq.gz \ --CPU 6 # 有基因组引导组装 Trinity --genome_guided_bam rnaseq_alignments.csorted.bam --max_memory 50G \ --genome_guided_max_intron 10000 --CPU 6 在运行中过程中,需要注意以下几点 质量控制(Quality control)。Trinity的--trimmomatic参数会调用Trimmomatic对数据进行过滤,这一步可以用其他软件完成。目前的RNA-seq质量也不需要过多的过滤。 Trinity中有一个"In silico Read Normalization",用于对read进行标准化,适用于超过300M的数据,默认开启,可以用--no_normalize_reads关闭。标准化的原因是,由于某些高表达基因会被检测到很多次,但是对于组装没有帮助,所以可以提前剔除。 如果基因组中基因密度大(比如说真菌),一些转录本可能会在UTR区域有重叠。那么为了尽可能降低转录本的错误融合,需要用到--jaccard_clip。对于植物和脊椎动物,就不需要考虑这一步。 运行结束后,最后的结果是trinity_out_dir的Trinity.fasta.Trinity将含有相同序列的转录本进行聚类,这种聚类可以被粗粗的被认为成一个基因的多个转录本。举个例子 >TRINITY_DN1000|c115_g5_i1 len=247 path=[31015:0-148 23018:149-246] AATCTTTTTTGGTATTGGCAGTACTGTGCTCTGGGTAGTGATTAGGGCAAAAGAAGACAC ACAATAAAGAACCAGGTGTTAGACGTCAGCAAGTCAAGGCCTTGGTTCTCAGCAGACAGA AGACAGCCCTTCTCAATCCTCATCCCTTCCCTGAACAGACATGTCTTCTGCAAGCTTCTC CAAGTCAGTTGTTCACAGGAACATCATCAGAATAAATTTGAAATTATGATTAGTATCTGA TAAAGCA "TRINITY_DN1000|c115" 是Trinity 聚类编号,“g5”是基因编号,“i1”是转录亚型编号 评估组装质量 有如下几种方法可以评估组装的质量 使用Bowtie/BWA将RNA-seq回贴到组装的转录组上,有80%以上的回帖率就行了。 用全长重构蛋白编码基因去搜索已知蛋白序列,见representation of full-length reconstructed protein-coding genes, 使用BUSCO根据保守同源基因进行评估 计算E90N50, 计算DETONATE得分 使用TransRate评估转录组组装

使用docker在CentOS7上搭建WordPress

自己课题组的网页挂了,老板一句话让我们自己写一个。结果被我拖了半年,结果丢给师弟,师弟决定用WordPress, 于是我就去给他配置服务器环境。 首先你得知道什么叫做WordPress, 它是一个基于PHP和MySQL的开源的博客管理工具,用于管理你的写作内容。由于它十分容易部署,而且有很多好看的主题可以供你选择,因此被许多人使用。 一般的安装方法是,你得有一个服务器,然后在服务器上按照PHP,MySQL, Apache/NGINX, 之后下载WordPress的安装包,进行编译安装。 只从有了docker,环境配置部分就得到了简化,并且你不用担心在准备环境的时候,要去调整原来的PHP版本,去修改MySQL的版本,还要专门折腾Apache或者NGINX。 后续的步骤都要求你对服务器有绝对的控制权,也就是说后续没有特别说明,我们都是以root用户运行命令 安装Docker Docker是目前服务器部署届的佼佼者,无论是部署网页,还是部署你的生信分析平台,只要你写好dockfile(一种描述部署规则的文件), 在任意的服务器上,安装好Docker,就可以构建出一个完全一样的运行环境。 我用的的是CentOS7, 版本信息如下 uname -r 3.10.0-862.el7.x86_64 然后用YUM工具进行安装, yum update -y yum install docker -y 由于国情,我们需要对配置一下docker的下载镜像,提高一下后续的加载速度。 使用vim编辑 /etc/docker/daemon.json, 增加如下内容。 "registry-mirrors": ["https://6xacs6l2.mirror.aliyuncs.com"] 启动我们的docker服务 systemctl start docker.service 默认情况下,只有root用户和docker组的用户才能访问Docker引擎的Unix socket。当然直接用root权限使用docker太过危险,建议新建一个docker用户,然后加入docker用户组 groupadd docker useradd -g docker docker 后续的操作用su docker切换到docker用户进行。 安装WordPress 首先让我们拉取WordPress的镜像 docker pull wordpress:latest WordPress的详细使用方法见https://hub.docker.com/_/wordpress/, 运行命令是 docker run --name some-wordpress --link some-mysql:mysql -d wordpress 先做一波参数说明: --name 容器的的名字 --link 和其他容器做连接 ·-d/--detach` 后台运行 但是直接这样运行是绝对会报错的,你需要再拉取一个MySQL容器, docker pull mysql:5.6 从MySQL镜像中运行单独的容器 docker run -d --privileged=true --name myMysql -v /data/mysql:/var/lib/mysql -e MYSQL_ROOT_PASSWORD=123456 -p 33306:3306 mysql:5.6 再来一波参数解释: -p: 端口映射,33306表示宿主,3306表示容器中的端口。 这里表示将宿主机的33306映射给镜像的3306. -e: 环境变量, 环境变量和具体的Docker容器制作时设置有关,这里表示设置镜像中MySQL的root 密码时123456 -v: 指定数据卷,也就是将我们MySQL容器的/var/lib/mysql映射到宿主机的/data/mysql --privileged=true: CentOS系统下的安全Selinux禁止了一些安全权限,导致MySQL容器在运行时会因为权限不足而报错,所以需要增加该选项 用docker ps -a 查看MySQL是否正常运行,出现错误的话,需要用docker stop 容器名停止运行,然后用docker rm 容器名删除容器,之后去掉-d选项重新运行排查错误。 当我们配置好MySQL后,下一步才是运行WordPress docker run -d --name mwp -e WORDPRESS_DB_HOST=mysql -e WORDPRESS_DB_PASSWORD=123456 -p 1080:80 --link myMysql:mysql wordpress 这里我们设置了两个环境变量,"WORDPRESS_DB_HOST"和"WORDPRESS_DB_PASSWORD", 但比较重要的有下面几个 "WORDPRESS_DB_HOST": 链接的docker的MySQL的IP地址和端口,一般设置成mysql表示用默认的设置 "WORDPRESS_DB_USER": 以什么用户使用MySQL,默认是root "WORDPRESS_DB_PASSWORD" 这设置MySQL的登陆用户密码,由于上一项是默认的root,所以这一项和之前的"MYSQL_ROOT_PASSWORD“要相同。 "WORDPRESS_DB_NAME": 数据库的表名,不需要修改,用默认的”wordpress"就行 之后在浏览器上用你服务器的IP,和映射出的端口号(我的是1080),就会得到配置界面 注意:尽管将容器的80端口映射给主机的1080,不需要用到root权限,但CentOS默认的防火墙禁止了大于1000后的所有端口,所以我们要开启这个端口 firewall-cmd --zone=public --add-port=8000/tcp --permanent firewall-cmd --reload #重载 后面就是你自己挑选主题的时刻了。 参考资料: https://www.cnblogs.com/qingyunzong/p/9011006.html https://blog.csdn.net/xiaoping0915/article/details/79515309 https://hub.docker.com/_/wordpress/

R和Rstudio的几个中文相关报错解决方案

果子老师做了一个非常详细的新手入门R语言的安装策略,叫做新手第1课,无敌无脑的R语言环境配置教程。基本上,你只要照着他的说的做,一字一句的阅读他的文档里的内容(注意,一定要一字一句),基本上R语言就能顺利用起来了。 不过在让我的一个师妹帮忙测试的时候,遇到了一些问题,这些问题无一例外的都和中文用户名有关。 这些问题基本上会在正版的win10出现,因为大家新买了一台笔记本都会预装了正版win10,然后启用的时候有些人会用中文名做用户名,所以就可能报错。 第一个问题:Fatal error: ERROR system error 5 系统属性-环境变量 我这里是字母作为用户名,所以没有问题,如果是中文那么就会报错。解决方案就是在随便一个盘中创建一个文件夹(不能有中文字符),然后把TEMP的路径改成该文件夹所在路径即可。 第二个问题:file.edit('~/.Rprofile') 没有打开一个新的文件。

使用purge_haplogs处理基因组杂合区域

基因组某些区域可能有着比较高的杂合度,这会导致基因组该区域的两个单倍型被分别组装成primary contig, 而不是一个为primary contig, 另一个是associated haplotig. 如果下游分析主要关注于单倍型,这就会导致一些问题。 那么有没有解决方案呢?其实也很好办,就是找到相似度很高的contig,将他们拆分。 purge_haplogs根据minimap2的比对结果,通过分析比对read的覆盖度决定谁去谁留。该工具适用于单倍型组装软件,例如 Canu, FALCON或 FALCON-Unzip primary contigs, 或者是分相后的二倍体组装(Falcon-Unzip primary contigs + haplotigs 。 purge_haplotigs依赖软件比较多,手动安装会很麻烦,但是他可以直接用bioconda装 conda create -n purge_haplotigs_env conda activate purge_haplotigs_env conda install purge_haplotigs 安装完成后需要一步测试 purge_haplotigs test 数据准备。 需要下载的数据集分为两个部分,一个是FALCON-Unzip后的primary contig 和 halplotigs. 另一个则是已经比完后的BAM文件 mkdir purge_haplotigs_tutorial cd purge_haplotigs_tutorial wget https://zenodo.org/record/841398/files/cns_h_ctg.fasta wget https://zenodo.org/record/841398/files/cns_p_ctg.aligned.sd.bam # 1.7G wget https://zenodo.org/record/841398/files/cns_p_ctg.aligned.sd.bam.bai wget https://zenodo.org/record/841398/files/cns_p_ctg.fasta wget https://zenodo.org/record/841398/files/cns_p_ctg.fasta.fai 当然我们不可能直接就拿到比对好的BAM文件,我们一般是有组装后的基因组以及用于组装的subread,假设这两个文件命名为, genome.fa 和 subreads.fasta.gz. minimap2 -ax map-pb genome.fa subreads.fasta.gz \ | samtools view -hF 256 - \ | samtools sort -@ 8 -m 1G -o aligned.bam -T tmp.ali 如果你有二代测序数据,也可以用BWA-MEM进行比对得到BAM文件。 第一步:使用purge_haplotigs readhist从BAM中统计read深度,绘制柱状图。 samtools mpileup -r "000005F|quiver" -f cns_p_ctg.fasta cns_p_ctg.aligned.sd.bam 也就是下图,你能明显的看到图中有两个峰,一个是单倍型的覆盖度,另一个二倍型的覆盖度, 第二步: 根据read-depth信息选择阈值。 purge_haplotigs contigcov -i cns_p_ctg.aligned.sd.bam.gencov -o coverage_stats.csv -l 20 -m 75 -h 190 这一步生成的文件是"coverage_stats.csv" 第三步:区分haplotigs. purge_haplotigs purge -g cns_p_ctg.fasta -c coverage_stats.csv -b cns_p_ctg.aligned.sd.bam -t 4 -a 60 这一步会得到如下文件 curated.artefacts.fasta:无用的contig,也就是没有足够覆盖度的contig. curated.fasta:新的单倍型组装 curated.haplotigs.fasta:从原本组装分出来的haplotigs curated.reassignments.tsv: 单倍型的分配信息 curated.contig_associations.log: 运行日志, 下面是其中一个记录,表示000004F_004和000004F_027是000004F_017的HAPLOTIG, 而000004F_017和000004F_013又是000004F,的HAPLOTIG。 000004F,PRIMARY -> 000004F_013,HAPLOTIG -> 000004F_017,HAPLOTIG -> 000004F_004,HAPLOTIG -> 000004F_027,HAPLOTIG 由于我们用的是单倍型组装primary contigs而不是二倍体组装的parimary + haplotigs, 因此我们需要将FALCON_Unzip的haplotgi合并到重新分配的haplotigs中,这样子我们依旧拥有二倍体组装结果 cat cns_h_ctg.fasta >> curated.haplotigs.fasta

官方提供了编译好的jar包,方便使用 wget https://github.com/broadinstitute/pilon/releases/download/v1.22/pilon-1.22.jar java -Xmx16G -jar pilon-1.22.jar 如果要顺利运行程序,要求JAVA > 1.7, 以及根据基因组大小而定的内存,一般而言是1M大小的基因对应1GB的内存。 Pilon有如下作用 对初步组装进行polish 寻找同一物种不同株系间的变异,包括结构变异检测 他以FASTA和BAM文件作为输入,根据比对结果对输入的参考基因组进行提高,包括 单碱基差异 小的插入缺失(indels) 较大的插入缺失或者block替换时间 填充参考序列中的N 找到局部的错误组装 最后它输出polish后的FASTA文件, 以及包含变异信息的VCF文件(可选) 推荐使用PCR-free建库测序得到的Illumina paired-end数据,这样子避免了PCR-duplication,有效数据更多,也不需要在分析过程中标记重复。 下面步骤,假设你的组装文件为draft.fa, 质量控制后的illumina双端测序数据分别为read_1.fq.gz和read_2.fq.gz 第一步:比对 bwa index -p index/draft draft.fa bwa mem -t 20 index/draft read_1.fq.gz read_2.fq.gz | samtools sort -@ 10 -O bam -o align.bam samtools index -@ 10 align.bam 第二步:标记重复(非PCR-free建库) sambamba markdup -t 10 align.bam align_markdup.bam 第三步:过滤高质量比对的read samtools view -@ 10 -q 30 align_markdup.bam > align_filter.bam samtools index -@ 10 align_filter.bam 第三步:使用Pilon进行polish MEMORY= #根据基因组大小而定 java -Xmx${MEMORY}G -jar pilon-1.22.jar --genome draft.fa --frags align_filer.bam \ --fix snps,indels \ --output pilon_polished --vcf &> pilon.log 一般Pilon迭代个2到3次就够了,所谓事不过三,过犹不及。 关于Pilon的一些参数说明: --frags表示输入的是1kb以内的paired-end文库,--jumps表示 大于1k以上的mate pair文库, --bam则是让软件自己猜测 -vcf: 输出一个vcf文件,包含每个碱基的信息 --fix: Pilon将会处理的内容,基本上选snps和indels就够了 --variant: 启发式的变异检测,等价于--vcf --fix all,breaks, 如果是polish不要使用该选项 minmq: 用于Pilon堆叠的read最低比对质量,默认是0。 阅读日志输出 这个日志文件是标准输出而不是标准错误输出,不过保险起见用&> 最开始,Pilon会输出他的版本信息(如下示例),以及将会对基因组做的调整, Pilon version 1.14 Sat Oct 31 14:30:00 2015 -0400 Genome: genome.fasta Fixing snps, indels 其中Fixing后面的含义为: "snps": 单碱基差异 "indels":小的indel的差异 "amb": 替换原有的N "gaps": 填充基因组的gap "local": 检测和修改错误组装 "all": 上述所有 "none": 不是上述的任何一种 接着Pilon会分染色体对BAM文件进行处理,根据BAM文件进行堆叠(pileup), 这个时候它会输出有效reads的深度,这里的有效reads包括未成对的read和正确成对的read。 Processing ctg1:1-5414473 frags align_mkdup.bam: coverage 19 Total Reads: 808985, Coverage: 19, minDepth: 5 从Pilon v1.4开始,Pilon还会输出基因组得到确认的碱基比例。 Confirmed 5403864 of 5414473 bases (99.80%) 后续是Pilon将会对原参考基因组做的一些调整的总体情况,如下表示纠正2个snp, 2个小的插入,4个缺失。 Corrected 2 snps; 0 ambiguous bases; corrected 2 small insertions totaling 12 bases, 4 small deletions totaling 6 bases 最后声明当前部分处理结束 Finished processing ctg1:1-5414473 如果,在--fix中选了gaps, 那么输出的内容还有如下内容。其中82048 -0 +276解释为在坐标82428处移除0个碱基,插入276个碱基。 # Attempting to fill gaps fix gap: scaffold00001:82428-93547 82428 -0 +276 ClosedGap https://github.com/broadinstitute/pilon/wiki/Standard-Output https://github.com/broadinstitute/pilon/wiki/Methods-of-Operation

都8102年了,还用fastq-dump,快换fasterq-dump吧

之前写过一篇文章Fastq-dump: 一个神奇的软件, 详细介绍了fastq-dump的用法。 虽然fastq-dump参数很多,而且一直被吐槽参数说明写的太差,但是如果真的要用起来其实也就是一行代码 fastq-dump --gzip --split-3 --defline-qual '+' --defline-seq '@$ac-$si/$ri' SRRXXXXX| SRRXXXX.sra # 加上--gzip后需要时间进行文件压缩 当然除了参数问题,还有一个让人诟病的地方就是他只能单个线程,所以速度特别的慢。尽管相对于下游分析要分析好几天而言,这点时间还能能等的。但是能快一点总是好的,所以在2018年的6月份,sra-tools更新了一个新的sra解压工具,fasterq-dump, a faster fastq-dump,它能利用临时文件和多线程加速从SRA文件提取FASTQ。 fasterq-dump的用法和fastq-dump一样,如下所示 fasterq-dump --split-3 ./SRR5318040 此外还有建立了GitHub Wiki提供使用教程,参见https://github.com/ncbi/sra-tools/wiki/HowTo:-fasterq-dump。 重点参数是-e|threads, 用于选择使用多少线程进行运行,默认是6个线程。 同时考虑到有些人容易着急,还提供了-p选项用于显示当前进度。 我用一个9G大小的SRA文件,分别以fastq-dump和fasterq-dump进行了测试。 time fastq-dump --split-3 -O test SRR5318040.sra # 558.76s user 41.36s system 101% cpu 9:51.82 total time fasterq-dump --split-3 ./SRR5318040 -e 20 -o SRR5318040 # 582.70s user 121.06s system 1130% cpu 1:02.25 total 从用户模式(user mode)来看, 两者的总CPU使用时间都差不多是560秒,从内核模式来看(Kernel Mode)来看,fasterq-dump花了更多时间在调用底层硬件上,例如分配内存地址。fastq-dump基本上稳定在一个线程,而fasterq-dump尽管指定了20个线程,但平均只用了11.5个线程吧。 对于我们而言,我们只要看最后的total部分,也就是实际花了多少时间。fastq-dump花了快10分钟,而fasterq-dump只需要1分钟,快了9倍多。 最后还有一点不足之处:输出的fastq的ID目前暂时没有选项可以调整,需要自己写个脚本解决。

读源码学C之阅读李恒的bioawk

目前尚没有能力直接去阅读htslib的源代码,看到bioawk的代码稍微简单点,因此准备先从这里下手,bioawk的项目地址为https://github.com/lh3/bioawk。 这次先阅读了"main.c"部分学习如何解析参数。 1-23: 版权信息部分 序列 指的是一组对象的集合,其中允许重复。序列分为有限序列和无限序列两种类型,我们通常用 表示序列中的第n个对象。 递归其实就是当前的序列依赖于之前的序列。最典型的案例就是兔子繁衍问题,假设一开始第一个月有1对兔子(A),到了第二个月这对兔子(A)性成熟,到了第三个月这对兔子(A)生出了一对兔子(B)。到了第四个月兔子(A)又生了一对兔子(D),之前刚的兔子(B)性成熟。到了第五个月,兔子(A)又生了兔子E,兔子(B)生出了兔子(F),而兔子(D)性成熟。 当然一图胜过千言万语,也就是下面这个情况 long int n = atoi(argv[1]); long int k = atoi(argv[2]); printf("month is %ld, %ld rabbit pairs per \n", n, k); long total = FIB(n,k); printf("The total number of rabbit is %ld \n", total); return 0; 注意,这个数字增长的很快,所以用了long int。 上面的方法在month很大的时候时候,这些迭代算法的由于会反复求解一些重复的子问题,导致求解速度就会非常慢,比如求解第100个月会有多少个兔子,会需要很久很久。而且很容易导致栈溢出,解决方案就是**动态规划**dynamic programming 动态规划在查找有很多重叠子问题的情况的最优解时有效。它将问题重新组合成子问题。为了避免多次解决这些子问题,它们的结果都逐渐被计算并被保存,从简单的问题直到整个问题都被解决。因此,动态规划保存递归时的结果,因而不会在解决同样的问题时花费时间。 代码如下: #include <stdio.h> #include <stdlib.h> int main(int argc, char *argv[]) if (argc != 2){ printf("USAGE: ./FIBv2 n k \n"); long int n = atoi(argv[1]); long int k = atoi(argv[2]); long int month[n-1]; int i = 0; for (i = 0 ; i < n; i++){ if ( i == 0 || i == 1){ month[i] = 1; } else{ month[i] = month[i-1] + k * month[i-2]; printf("After %ld month, there are %ld rabbits\n", n, month[n-1]); return 0; 这个代码在求解第100个月有多少兔子的时候,基本上是秒答。

整合了串联重复序列和遍在重复序列的屏蔽(之前没有这一步) 以GFA格式存放graph文件,后续可以用Bandage进行可视化 通过算法和性能优化,提高了Associate Contigs的准确性 分析流程的性能优化 软件安装和数据准备 Falcon终于拥抱了bioconda, 这也就意味着我们再也不需要用到他们原本笨拙的安装脚本,浪费时间在安装软件上。 conda create -n pb-assembly pb-assembly source activate pb-assembly conda create -p ~/opt/biosoft/pb-assembly pb-assembly source activate ~/opt/biosoft/pb-assembly 这里使用https://pb-falcon.readthedocs.io/en/latest/tutorial.html上的所用的E. coli数据集 wget https://downloads.pacbcloud.com/public/data/git-sym/ecoli.m140913_050931_42139_c100713652400000001823152404301535_s1_p0.subreads.tar.gz tar -xvzf ecoli.m140913_050931_42139_c100713652400000001823152404301535_s1_p0.subreads.tar.gz 解压缩后的文件夹里有三个300M的Fasta文件, 将他们的实际路径记录到input.fofn中 ecoli.1.fasta ecoli.2.fasta ecoli.3.fasta 准备配置文件 为了进行组装,需要准备一个配置文件。我的配置文件为fc_run.cfg,内容如下。你们可以先预览一下,后面看我的解释说明。 #### Input [General] input_fofn=input.fofn input_type=raw pa_DBdust_option= pa_fasta_filter_option=pass target=assembly skip_checks=False LA4Falcon_preload=false #### Data Partitioning pa_DBsplit_option=-x500 -s50 ovlp_DBsplit_option=-x500 -s50 #### Repeat Masking pa_HPCTANmask_option= pa_REPmask_code=1,100;2,80;3,60 ####Pre-assembly genome_size=0 seed_coverage=20 length_cutoff=3000 pa_HPCdaligner_option=-v -B128 -M24 pa_daligner_option=-e.7 -l1000 -k18 -h80 -w8 -s100 falcon_sense_option=--output-multi --min-idt 0.70 --min-cov 2 --max-n-read 800 falcon_sense_greedy=False ####Pread overlapping ovlp_daligner_option=-e.96 -l2000 -k24 -h1024 -w6 -s100 ovlp_HPCdaligner_option=-v -B128 -M24 ####Final Assembly overlap_filtering_setting=--max-diff 100 --max-cov 300 --min-cov 2 fc_ovlp_to_graph_option= length_cutoff_pr=2000 [job.defaults] job_type=local pwatcher_type=blocking JOB_QUEUE=default MB=32768 NPROC=6 njobs=32 submit = /bin/bash -c "${JOB_SCRIPT}" > "${JOB_STDOUT}" 2> "${JOB_STDERR}" [job.step.da] NPROC=4 MB=32768 njobs=32 [job.step.la] NPROC=4 MB=32768 njobs=32 [job.step.cns] NPROC=8 MB=65536 njobs=5 [job.step.pla] NPROC=4 MB=32768 njobs=4 [job.step.asm] NPROC=24 MB=196608 njobs=1 根据注释信息,文件分为"input", "Data partitioning", "Repeat Masking", "Pre-assembly", "Pread overlapping", "Final Assembly", 以及最后的任务调度部分,让我们分别看下这里面的内容 输入(Input) #### Input [General] input_fofn=input.fofn input_type=raw pa_DBdust_option= pa_fasta_filter_option=pass target=assembly skip_checks=False LA4Falcon_preload=false 输入这里参数比较简单,基本不需要做任何改动,除了 pa_fasta_filter_options,用于处理一个ZMW(测序翻译孔)有多条subread时,到底选择哪一条的问题。 "pass": 不做过滤,全部要。 "streamed-median": 表示选择大于中位数的subread "streamed-internal-median": 当一个ZMW里的subread低于3条时选择最长,多于单条则选择大于中位数的subread 0.01版本pb-assembly的pa_DBdust_option有一个bug,也就是里面的参数不会传递给DBdust, DBdust是对read进行soft-masking,一般都用默认参数,因此这个bug问题不大。 数据分配(Data Partitioning) pa_DBsplit_option=-x500 -s50 ovlp_DBsplit_option=-x500 -s50 这部分的设置会将参数传递给DBsplit,将数据进行拆分多个block,后续的并行计算都基于block。-s 50表示每个block大小为50M。 这适用于基因组比较小的物种,如果是大基因组则应该设置为-s 200或者-s 400 重复屏蔽(Repeat Masking) #### Repeat Masking pa_HPCTANmask_option= pa_REPmask_code=1,100;2,80;3,60 屏蔽重复序列可以在不损失组装准确性的同时,提高后续组装的overlap/daligner步骤10~20倍速度,见Detecting and Masking Repeats. pa_HPCTANmask_option的参数会传给串联重复步骤的HPCTANmask, 而pa_REPmask_code很复杂,它分为三次迭代,因此这里1:100;2,80;3,60 就表示第一次迭代检测每个block中出现超过100次的序列,第二次迭代将2个block合并一起检测超过80次的序列,第三次将3个block进行合并检测超过60次的序列。 https://github.com/PacificBiosciences/pbbioconda/issues/20 预组装(纠错)pre-assembly ####Pre-assembly genome_size=0 seed_coverage=20 length_cutoff=3000 pa_HPCdaligner_option=-v -B128 -M24 pa_daligner_option=-e.7 -l1000 -k18 -h80 -w8 -s100 falcon_sense_option=--output-multi --min-idt 0.70 --min-cov 2 --max-n-read 800 falcon_sense_greedy=False 当length_cutoff=-1时,设置genome_size和seed_coverage会自动计算要过滤的序列。否则是过滤低于一定长度的read。 pa_HPCdalinger_option参数不需要调整,-M表示每个进程的内存为24G,一般200M的block对应16G。 pa_daligner_option的参数比较重要: -e:错误率,低质量序列设置为0.70,高质量设置为0.80。 值越高避免单倍型的坍缩 -l: 最低overlap的长度,文库比较短时为1000, 文库比较长为5000. -k: 低质量数据为14,高质量数据为18 -h: 表示完全match的k-mer所覆盖的碱基数。和-l, -e有关,越大越严格。预组装时最大也不要超过最低overlap长度的1/4. 最低就设置为80 纠错后相互比对 ####Pread overlapping ovlp_daligner_option=-e.96 -l2000 -k24 -h1024 -w6 -s100 ovlp_HPCdaligner_option=-v -B128 -M24 和上面的参数类似,但是-e的范围调整为0.93-0.96,-l范围调整为1800-6000, -k调整为 18-24 overlap_filtering_setting=--max-diff 100 --max-cov 300 --min-cov 2 fc_ovlp_to_graph_option= length_cutoff_pr=2000 这里的参数就可以随便调整了,因为这一步速度很快。 例如length_cutoff_pr就可以从2000,提高到15000. 最后还有一部分是任务投递系统,如果是单节点运行,需要注意设置 njobs,这是同时投递的任务数。假如你将[job.step.cns]按照如下的方式设置,那么同时会出现 个任务,如果你的内存只有128G,运行一段时间后你的所有内存就会被耗尽,那么基本上你就只能重启服务器了。 [job.step.cns] NPROC=8 MB=65536 njobs=50 用上述的配置文件,以fc_run fc_run.cfg运行后,最后的2-asm-falcon/p_ctg.fa的序列数有4条,最长为4,685,024, 之后我将length_cutoff_pr调整为15k,2-asm-falcon/p_ctg.fa序列只有一条,长度为4,638,411 下载asm.gfa到本地,用Bandage可视化,可以发现组装效果不错。 注意:这一步必须在后续步骤之前运行。 通常,我们需要准备一个物种的基因组fasta文件,当然RNA和protein都是没有问题。通过prepare-refseqs.pl格式化生成的track,这为后续所有文件提供一个坐标,一直放大后参考序列的碱基也会显示出来。生成的track 会为后续所有文件提供一个坐标,一直放大后参考序列的碱基也会显示出来。 主要用到工具是prepare-refseqs.pl,他的用法很多,如下: prepare-refseqs.pl --gff <GFF file> [options] # OR: prepare-refseqs.pl --fasta <file1> --fasta <file2> [options] # OR: prepare-refseqs.pl --indexed_fasta <file> [options] # OR: prepare-refseqs.pl --twobit <file> [options] # OR: prepare-refseqs.pl --conf <JBrowse config file> [options] # OR: prepare-refseqs.pl --sizes <sizes file> [options] 更多内容用prepare-refseqs.pl -h查看。 # 以下操作在jbrowse家目录,序列文件根据实际情况修改 bin/prepare-refseqs.pl --fasta ~/lyrata/Sequence/Alyrata_384_v1.fa 这就会在当前生成data文件夹,直接访问I地址所看到的序列就来源于该文件夹。 准备特征序列 特征序列一般以"gff|gbk|bed"格式存放,用于注明序列的信息。所需工具为flatfile-to-json.pl bin/flatfile-to-json.pl --gff ~/lyrata/Annotation/Alyrata_384_v2.1.gene.gff3 --trackType CanvasFeatures --trackLabel lyrata 结果是在当前目录下生成data,data里包括序列track配置文件. 同样可以用--out参数输出到指定文件夹。 bin/flatfile-to-json.pl --gff ~/Athalina/TAIR10/TAIR10_GFF3_genes.gtf --tracklabel gene --out ./Athaliana/ 除了在JBrowse上通过具体位置定位外,我们还可以0JBrowse上通过基因名快速定位到目标区间,只需要在上两步的基础上运行下面程序即可。 bin/generate-names.pl 推荐阅读: http://jbrowse.org/docs/prepare-refseqs.pl.html http://jbrowse.org/docs/generate-names.pl.html

JBrowse使用说明:如何安装JBrowse

JBrowse is a fast, scalable genome browser built completely with JavaScript and HTML5. It can run on your desktop, or be embedded in your website. 如果你想要使用JBrowse,一定要有管理员权限,否则建议使用IGV。 你的服务器必须安装有libpng,zlib,libgd,make,C compiler, C++ compiler。 对不同的操作系统有不同的安装方法,如果你有管理员权限,那么用系统自带的管理工具就行。 # Ubuntu/Debian sudo apt-get install build-essential libpng-dev zlib1g-dev libgd2-xpm-dev # Red Hat/Fedora/CentOS sudo yum groupinstall "Development Tools" sudo yum install libpng-devel gd-devel zlib-devel perl-ExtUtils-MakeMaker # Mac OS X (homebrew) ## 需要在APP Store 安装 xcode 或在命令行输入gcc提示安装 xcode-select --install brew install libpng libgd zlib 注:如果你没有管理员权限,就请自行编译并添加到环境变量中。你会发现这非常地麻烦。 安装jbrowse jbrowse有两种安装形式: 系统级: 面向更多用户,常常部署到网页服务器 普通用户级:满足内网用户的一般需求 系统级安装(最高权限) 系统级安装(需要管理员权限) # make a directory that this user can write to # for ubuntu/Debian that is /var/www/; for centos is /var/www/html sudo mkdir /var/www/jbrowse; sudo chown `whoami` /var/www/jbrowse; # cd into it cd /var/www/jbrowse; # fetch a JBrowse release zip file curl -O http://jbrowse.org/releases/JBrowse-1.12.3/JBrowse-1.12.3.zip # unzip it and cd into it unzip JBrowse-1.12.3.zip cd JBrowse-1.12.3 # run setup.sh, quick start with example data ./setup.sh 系统级别的JBrowse安装到/var/www,然后可以通过http://机器IP地址/jbrowse/JBrowse-1.12.3/docs/tutorial/index.html判断安装是否成功。 但是,如果你的服务器上没有安装"apache"或"nginx"等web服务器软件时,上述操作无法成功访问网页。 以下句子你可能已经看不懂了: 尽管软件已经部署在你的服务器上,但是仅仅是你本人才能查看或使用软件,这对外人不可见。就像你家买了一台老罗的畅呼吸,你不和别人说,别人是不知道的。如果他需要你家这台畅呼吸解决他家的空气问题,”唰的一下就没了“,就需要专门铺设一个管道连接你们两家,这样子你家的畅呼吸就把他家的空气也净化了一下。 因此,你需要安装一个web服务器软件,然后设置规则,让别人能够顺利访问到你的本地资源。这里安装的是nginx. # 安装nginx,基于ubuntu sudo apt-get install nginx # 检查是否安装成功 sudo nginx -t nginx检查 以上网页存放在/var/www/html/index.nginx-debian.html,由/etc/nginx/sites-enabled/default定向所在目录。其中/etc/nginx/sites-enabled/default是/etc/nginx/sites-available/default的软连接。这是nginx安装完成之后默认设置,提供一个案例。如果我们需要配置nginx,使其找到Jbrowse所在目录,就需要在/etc/nginx/sites-available/创建配置文件,然后链接到/etc/nginx/sites-enabled/. sudo vi /etc/nginx/sites-avaiable/jbrowse # 内容如下 server { listen 8080; listen [::]:8080; root /var/www/jbrowse/JBrowse-1.12.3; index index.html # 链接到/etc/nginx/sites-enabled/ sudo ln -n /etc/nginx/sites-available/jbrowse /etc/nginx/sites-enabled/jbrowse ps -ef | grep nginx sudo kill -QUIT nginx 主进程号 sudo nginx 然后你就发现直接可以通过IP:8080就能访问JBrowse,而不是要求那样的一长串地址。 普通用户级(不需要管理员权限):和系统级的差异,文件存放在用户目录下。 # 软件原始文件,我一般存放在家目录的src文件夹下 cd ~/src &curl -O http://jbrowse.org/releases/JBrowse-1.12.3/JBrowse-1.12.3.zip # 解压缩,并移动到我软件文件夹下 unzip JBrowse-1.12.3.zip && mv JBrowse-1.12.3 ~/biosoft cd ~/biosoft/JBrowse-1.12.3 && ./setup.sh 普通用户:nginx和apche是针对大型访问站点,如果你只是实验室内容访问,用Python做web服务器就行了。 cd ~/biosoft/JBrowse-1.12.3 python -m SimpleHTTPServer 5000 # 长时间运行 # nohup python -m SimpleHTTPServer 5000 & 是的就是如此简单,你就可以访问"IP:5000"访问页面。但是距离真实情况下的使用,还需要进行后续的配置 其实这两者在安装上没有多大区别,安装到/var/www下是因为文档编写者默认你的服务器已经有web服务器应用,并且那个应用托管了/var/www,而/var/www需要一定的权限。 但是实际上,如果你修改web服务器应用的配置,你可以将JBrowse安装在任意地方。甚至,对于少量的请求,可以直接用Python自带的web服务器让你的JBrowse被外部人员访问到。 Nginx从听说到学会 Nginx wiki 玩Python之HTTP代理 一起写一个 Web 服务器(1) 一起写一个 Web 服务器(2)

JBrowse使用说明:如何配置track分面搜索选择器

但JBrowse要展示数据达到一定量级之后,如何方便管理这些track就成了一个问题,JBrowse支持三种track选择器: JBrowse/View/TrackList/Simple:已经作古的选择方法 JBrowse/View/TrackList/Hierarchical:默认通过勾选的方法选择需要展示的track JBrowse/View/TrackList/Faceted:本文要介绍的分面搜索选择器(Track faceted track selector),如下图所示 然后新增一个bigwig track,之后以这个track为例进行讲解 bin/add-bw-track.pl --label 2cell --bw_url data/test.bw 对data/trackList.json进行修改,主要是增加“metadata"行 这里的jq需要去https://stedolan.github.io/jq下载安装。 cat trackList.json| jq -r '.tracks[] | [.label,.key] | @csv' > trackMetadata.csv 生成的csv文件可以下载到本地,用Excel进行调整,调成完成后覆盖原文件就行了,例如下面这个情况

Rosalind: DNA核苷酸计数和DNA翻译成RNA

DNA核苷酸计数 问题描述: 给定一行核苷酸序列,长度最长为1000 nt, 返回其中'A', 'T', 'C', 'G'出现的次数 C代码如下: #include <stdio.h> #include <string.h> #include <ctype.h> #include "dbg.h" #define MAX_SIZE 1001 int main( int argc, char *argv[]) check(argc == 2, "Need an argument"); FILE *filein ; if( ( filein = fopen(argv[1], "r")) == NULL){ printf("unable to open %s \n", argv[1]); return 1; char nucl[MAX_SIZE]; char *in = NULL; in = fgets(nucl, MAX_SIZE - 1, filein); check(in != NULL, "Failed to read file"); // count int count =0; char letter; int A_NUM = 0, T_NUM = 0, C_NUM = 0 , G_NUM = 0; for ( count = 0; nucl[count] != '\0' && nucl[count] != '\n'; count++){ letter = toupper(nucl[count]); switch(letter){ case 'A': A_NUM ++; break; case 'T': T_NUM ++; break; case 'C': C_NUM ++; break; case 'G': G_NUM ++; break; default: printf("%d: %c is not a valid base\n", count, letter); printf("A\tC\tG\tT\n"); printf("%d\t%d\t%d\t%d\t\n",A_NUM, C_NUM, G_NUM, T_NUM); fclose(filein); return 0; error: return -1; DNA转录成RNA 问题描述: RNA由A, C, G, U组成,其中U是对应DNA编码链上的T。给定一段编码链上的序列,将其翻译成RNA序列。 C代码如下: #include <stdio.h> #include <ctype.h> #include "dbg.h" #define LINE_MAX 80 int main(int argc, char *argv[]) check(argc == 2, "USAGES: ./RNA file"); char *filename = argv[1]; FILE *in = fopen(filename, "r"); check(in != NULL, "unable to open file"); char line[LINE_MAX]; int i = 0; while ( (fgets(line, LINE_MAX - 1, in)) != NULL){ for (i = 0; line[i] != '\0'; i++){ int base = toupper(line[i]); if ( base == 'T'){ printf("%c", 'U'); } else{ printf("%c", base); printf("\n"); return 0; error: return -1;

题图来源: https://pexels.com 只要开始学习,就会出现疑问。即便一个作者认为自己讲的多么具体,但是由于“知识的诅咒”存在,必然有一些点是他认为理所当然,但是读者/听众却从未听过的现象。譬如,之前在生信媛举办的一期互课活动中(互课:你的知识就是你的入场券),当时的主题是ATAC-seq数据分析,一些做遗传学的小伙伴就不知道什么叫做染色体开放区。 当我们遇到问题的时候,我们就会寻求解答。如果是学生时代,我们一般会选择向老师发问。如果有一个学霸是你隔壁,那你就会找他提问。现在搜索引擎极度发达,我们会选择上网检索,比如去知乎上提问,“有一个漂亮的女朋友是什么样的体验?”,“长得帅是一种什么样的体验?”。 这貌似提问是一件非常简单的小事,好像谁都会,但是其实结果千差万别。我经常在各种群里面见到的一种提问形式是,有人做XX吗?有一个问题想问下? image 看到这种问题,我内心的想法其实是“没有,滚”。这种问题仿佛是,多年没有联系的小学同学,突然有一天问你“在吗?”,你小心翼翼的点开他的朋友圈,发现半年前他领证了,经过强大的逻辑推理,你猜测他要办喜酒了,所以不要回复,不要回复,不要回复,你只要不回复他就不知道你在不在。 你永远不知道你回复“是”之后是什么结果。无厘头一点,如果有人问,“群里的大神,谁会折纸飞机吗?”,你回答“我会呀”,然后他接着一句,“我最近想造一架飞机,你能帮我吗?”。你是不是有一种想把那个提问者打爆的冲动。 除了上面“有人会做XX,我有一个问题想问下”体以外,我还见过直接丢出截图体,直接复制报错信息自己以为有用****的一行体,从我的经验看,这些问题的解决概率,就和你在路边拿个碗别人往你碗里丢钱的概率一样的,依赖于别人是否有心情搭理你。我也建议,看到这类提问,大家都不要回答,直到这些人学会了恰当的提问方式。 那你可能会问了,什么是恰当的提问方式呢?判断一个提问是否恰当的一个金标准就是,你是否在这个问题上花了足够多时间,让回答者看到你的诚意。以此延伸出的具体细节可以参考 提问的智慧: https://github.com/ryanhanwu/How-To-Ask-Questions-The-Smart-Way/blob/master/README-zh_CN.md 其中有三条非常重要: 只要是搜索引擎能回答的就别问别人。比如说你问别人如何将SAM文件转成BAM格式。这个问题你用SAM和BAM作为关键字出来一堆回答好吗? 你要知道自己提问的目的是什么?你会不会在作弊传递字条的时候,在字条写“在吗?”, 你肯定是要问,第1题第2问答案是什么吧!搞清楚自己的问题是什么,把问题问的清楚点,节约双方的时间。 要方便别人回答。所有人的时间都是有价值的,我不希望自己在解决别人的问题时候,还要推测出提问的人意图。你忘了你中学时被出题人的意图支配的恐惧了吗? 作为一个经常写点分析的教程的人,我也每天遇到问题,大部分问题我都通过百度、bing和谷歌解决了。还有一些问题,我暂时搞不定,我思索了半天,依旧有一些困惑,我会选择发邮件问下作者。举一个最近的一个案例吧, 我在摸索别人一句话带过的peak注释步骤时,有一步我存在疑惑,“terminal region as -1,000 to +1,000 bp from the poly(A) site”,我一直以来都用的是ChIPseeker和HOMER对peak进行注释,结果里没有出现的terminal region。为了避免我自己对GTF注释文件的不熟悉,我还特意去对GTF的第三列做了一下分析, 发现没有Poly(A) site的定义。虽然这里面也没有TSS,但是我们知道第一个外显子的位置就是TSS。带着疑问,我就发了一份邮件给对方。 由于他不是通讯作者,我并没有直接发邮件给通讯作者要他的联系方式,而是用通过各种途径去搜索(先找到了他的谷歌学术页面,发现他是通过学校邮箱验证的,于是用学校加他的中文拼音,最后找到了邮件地址),这就是能自己搜到的就别麻烦别人了。 我在自己薄弱的地方尽可能做了一些调查,在提问中做出了思考,并且将关键字粗体标明方便回答者找到,最后我也得到了我想要的答案。 其实除了这个问题外,我还有很多问题是原本要问但是没有问的。比如作者的peak calling这部分分析,我刚开始的时候得到peak size的中位数比他的大了200bp以上,和他原文的结论有出入。我本来是想发邮件问作者具体的脚本是什么?但是那个时候的我只跑了一种参数,其实还有很多可能性我没有运行,这就说明我还没有做足功课。于是我就继续尝试了多种策略,最后发现需要链特异性的peak calling需要在MACS2上加上一个—nomodel才行。 综上,让自己觉得问心无愧的提问才是一个合格的提问。

R-loop数据分析之R-ChIP(peak calling)

Peak Calling 关于MACS2的使用方法, 我写了如何使用MACS进行peak calling详细地介绍了它的参数,在用MACS2之前尽量去阅读下。 尽管 文章说他按照默认参数分别找narrow peak 和 broad peak, 即"We used MACS2 with default settings to call narrow (or broad when necessary) R-loop peaks",但是 MACS2的默认参数一开始会根据reads在正反链的分布建立双峰模型,确定偏移模型(shifting model),也就是它会认为示例CHTF8和CIRH1A位于正反链的信号视作一个IP结果的信号,因此用默认参数绝对有问题,所以 必须要增加--nomodel取消建模,且增加--shift 0 --extsize 150按照实验的条带长度对片段进行延长。 cd analysis # HKE293 D210N macs2 callpeak --nomodel --shift 0 --extsize 150 -g hs --keep-dup all -f BAM -t 2-read-align/HKE293-D210N-V5ChIP-Rep*.flt.bam -c 2-read-align/HKE293-D210N-Input-Rep*.flt.bam --outdir 5-peak-calling/narrow -n HKE293-D210 2> log/HKE293-D210.macs2.narrow.log & macs2 callpeak --nomodel --shift 0 --extsize 150 -g hs --keep-dup all --broad -f BAM -t 2-read-align/HKE293-D210N-V5ChIP-Rep*.flt.bam -c 2-read-align/HKE293-D210N-Input-Rep*.flt.bam --outdir 5-peak-calling/broad -n HKE293-D210 2> log/HKE293-D210.macs2.broad.log & # HEK293 WKKD macs2 callpeak --nomodel --shift 0 --extsize 150 -g hs --keep-dup all -f BAM -t 2-read-align/HKE293-WKKD-V5ChIP.flt.bam -c 2-read-align/HKE293-WKKD-Input.flt.bam --outdir 5-peak-calling/narrow/ -n HKE293-WKKD 2> log/HKE293-WKKD.macs2.narrow.log & macs2 callpeak --nomodel --shift 0 --extsize 150 -g hs --keep-dup all --broad -f BAM -t 2-read-align/HKE293-WKKD-V5ChIP.flt.bam -c 2-read-align/HKE293-WKKD-Input.flt.bam --outdir 5-peak-calling/broad/ -n HKE293-WKKD 2> log/HKE293-WKKD.macs2.broad.log & 可以将建模和不建模的narrowPeak结果在IGV进行比较,你会发现之前的CHTF8和CIRH1A上的两个信号峰在默认参数下就被认为是成一个。 文章中对找到peak进行了一次筛选,标准是大于5倍富集和q-value 小于或等于0.001(broad peak则是0.0001),最后文章写着在D210N有12,906peak,然后剔除WKKD里的peak,还有12,507个。 我找到的HKE293-D210的原始narrowPeak数为17,558, 按照作者的标准筛选后只剩下9,552,。对于HKE2930-D210的原始broadPeak数为42,912,过滤之后只剩2,998了。隐隐觉得peak有点太少了,于是我的标准是4倍变化。 负对照的WKKD无论是narrowPeak还是broadPeak都是0 cd analysis/5-peak-calling/ # narrow peak awk -v fc=$fc '$7 >= fc && $9 >=3' narrow/HKE293-D210_peaks.narrowPeak > HKE293-D210_flt.narrowPeak # broad peak awk -v fc=$fc '$7 >= fc && $9 >=4' broad/HKE293-D210_peaks.broadPeak > HKE293-D210_flt.broadPeak 之后可以用过滤后的D210N的narrowPeak和broadPeak的peak长度进行描述性统计分析,然后用箱线图展示其大小分布。 library(data.table) library(ggplot2) narrowPeak <- fread(file="HKE293-D210_flt.narrowPeak", sep="\t", header = F) broadPeak <- fread(file="HKE293-D210_flt.broadPeak", sep="\t", header = F) peak_size <- log10(c(narrowPeak$V3 - narrowPeak$V2, broadPeak$V3 - broadPeak$V2 )) peak_from <- factor(rep(c('Narrow','Broad'), times=c(nrow(narrowPeak),nrow(broadPeak)) ), levels=c("Narrow","Broad")) peak_df <- data.frame(size=peak_size, from=peak_from) my_clear_theme = theme_bw() + theme(panel.grid.major = element_line(colour="NA"), panel.grid.minor = element_line(colour="NA"), panel.background = element_rect(fill="NA"), panel.border = element_rect(colour="black", fill=NA), legend.background = element_blank()) ggplot(peak_df,aes(x=from,y=size,col=from)) + geom_boxplot(notch = T,outlier.size = .5) + scale_y_continuous(breaks=c(log10(100),log10(300),log10(400),log10(450),log10(500),log10(1000),log10(2000)), labels = c(100,300,400,450,500,1000,2000)) + coord_cartesian(ylim=c(2,3.5)) + labs(col="Peak Size") + my_clear_theme + xlab("") + ylab("Peak Size (bp)") ggsave("Fig2A_boxplot.pdf", width=4,height=6,units = "in") boxplot of peak width 原文说自己的narrow peak 长度的中位数是199bp, broad peak的长度中位数是318 bp,和电镜观察的150–500 bp一致。我过滤后的peak的中位数是208,broad peak的中位数是581。 如果你无脑用MACS2的参数话,就会得到右边的图,peak的长度明显都偏大了。 还可以对Broad peak 和Narrow peak进行比较,看看有多少共同peak和特异性peak。 #!/bin/bash peak1=${1?first peak} peak2=${2?second peak} outdir=${3?output dir} mkdir -p ${outdir} bedtools intersect -a $peak1 -b $peak2 > ${outdir}/$(basename ${peak1%%.*}).common.peak common=$(wc -l ${outdir}/$(basename ${peak1%%.*}).common.peak) echo "$common" bedtools subtract -A -a $peak1 -b $peak2 > ${outdir}/$(basename ${peak1}).only.bed left=$(wc -l ${outdir}/$(basename ${peak1}).only.bed) echo "$left" bedtools subtract -A -b $peak1 -a $peak2 > ${outdir}/$(basename ${peak2}).only.bed right=$(wc -l ${outdir}/$(basename ${peak2}).only.bed) echo "$right" 运行:bash scripts/09.peak_comparison.sh analysis/5-peak-calling/HKE293-D210_flt.narrowPeak analysis/5-peak-calling/HKE293-D210_flt.broadPeak analysis/5-peak-calling/peak_compare > results/narrow_vs_broad.txt 然后得到的数值就可以丢到R语言画图了,这个结果和Figure S2A一致。 library(VennDiagram) grid.newpage() venn.plot <- VennDiagram::draw.pairwise.venn(area1=6401 + 7165, area2=286 + 7165, cross.area = 7165, category = c("Narrow specific","Broad specific"), scaled = F, fill=c("#c7d0e7","#f4babf"), lty='blank', cat.pos = c(180,0),#360度划分 rotation.degree = 90 #整体旋转90 grid.draw(venn.plot) 我们可以对这三类peak对应区间内的信号进行一下比较。统计信号强度的工具是bigwigSummary,来自于ucscGenomeBrowser工具集。 脚本为scripts/10.peak_signal_summary.sh #!/bin/bash bed=${1?bed file} fwd=${2?forwad bigwig} rvs=${3?reverse bigwig} exec 0< $bed while read region chr=$( echo $region | cut -d ' ' -f 1) sta=$( echo $region | cut -d ' ' -f 2) end=$( echo $region | cut -d ' ' -f 3) (bigWigSummary $fwd $chr $sta $end 1 ; bigWigSummary $rvs $chr $sta $end 1 )| awk -v out=0 '{out=out+$1} END{print out/2}' 分别统计三类peak的信号的信号 bash scripts/10.peak_signal_summary.sh analysis/5-peak-calling/peak_compare/HKE293-D210_flt.broadPeak.only.bed analysis/4-normliazed-bw/HKE293-D210N-V5ChIP_fwd.bw analysis/4-normliazed-bw/HKE293-D210N-V5ChIP_rvs.bw > results/HKE293-D210N-V5ChIP.broadSignal bash scripts/10.peak_signal_summary.sh analysis/5-peak-calling/peak_compare/HKE293-D210_flt.narrowPeak.only.bed analysis/4-normliazed-bw/HKE293-D210N-V5ChIP_fwd.bw analysis/4-normliazed-bw/HKE293-D210N-V5ChIP_rvs.bw > results/HKE293-D210N-V5ChIP.narrowSignal bash scripts/10.peak_signal_summary.sh analysis/5-peak-calling/peak_compare/HKE293-D210_flt.common.peak analysis/4-normliazed-bw/HKE293-D210N-V5ChIP_fwd.bw analysis/4-normliazed-bw/HKE293-D210N-V5ChIP_rvs.bw > results/HKE293-D210N-V5ChIP.commonSignal 在R语言中绘图展示 # peak signal commonSignal <- fread("./HKE293-D210N-V5ChIP.commonSignal") narrowSignal <- fread("./HKE293-D210N-V5ChIP.narrowSignal") broadSignal <- fread("./HKE293-D210N-V5ChIP.broadSignal") data <- data.frame(signal=c(commonSignal$V1, narrowSignal$V1, broadSignal$V1), type=factor(rep(c("common","narrow","broad"), times=c(nrow(commonSignal),nrow(narrowSignal), nrow(broadSignal)))) ) wilcox.test(data[data$type=="common",1],data[data$type=="broad",1]) wilcox.test(data[data$type=="common",1],data[data$type=="narrow",1]) p <- ggplot(data,aes(x=type,y=signal,col=type)) + geom_boxplot(outlier.colour = "NA") + coord_cartesian(ylim=c(0,3.5)) + theme(panel.grid.major = element_line(colour="NA"), panel.grid.minor = element_line(colour="NA"), panel.background = element_rect(fill="NA"), panel.border = element_rect(colour="black", fill=NA)) + theme(axis.text.x = element_blank()) + xlab("") + ylab("R-loop Signal") + labs(col="") p1 <- p + annotate("segment", x=1,xend=1,y=0.75,yend=2,size=0.5) + annotate("segment", x=2,xend=2,y=1.25,yend=2,size=0.5) + annotate("segment", x=1,xend=2,y=2,yend=2,size=0.5) + annotate("text", x= 1.5, y=2.1, label="***") p2 <- p1 + annotate("segment", x=3,xend=3,y=0.65,yend=3,size=0.5) + annotate("segment", x=2,xend=2,y=2.2,yend=3,size=0.5) + annotate("segment", x=2,xend=3,y=3,yend=3,size=0.5) + annotate("text", x= 2.5, y=3.2, label="***")

R-loop数据分析之R-ChIP(样本间BAM比较和可视化)

样本间相关性评估 上一步得到各个样本的BAM文件之后,就可以在全基因组范围上看看这几个样本之间是否有差异。也就是先将基因组分成N个区间,然后用统计每个区间上比对上的read数。 脚本scripts/07.genome_bin_read_coverage.sh如下 #!/bin/bash set -e set -u set -o pipefail samples=${1?missing sample file} chromsize=${2:-index/hg19.chrom.sizes} size=${3:-3000} ALIGN_DIR="analysis/2-read-align" COV_DIR="analysis/3-genome-coverage" mkdir -p ${COV_DIR} exec 0< $samples while read sample bedtools makewindows -g $chromsize -w $size | \ bedtools intersect -b ${ALIGN_DIR}/${sample}.flt.bam -a - -c -bed > ${COV_DIR}/${sample}.ReadsCoverage 准备一个输入文件存放待处理样本的前缀,然后运行脚本bash scripts/07.genome_bin_read_coverage.sh samples_rep.txt HKE293-D210N-Input-Rep1 HKE293-D210N-Input-Rep2 HKE293-D210N-Input-Rep3 HKE293-D210N-V5ChIP-Rep1 HKE293-D210N-V5ChIP-Rep2 HKE293-D210N-V5ChIP-Rep3 最后将得到的文件导入到R语言中进行作图,使用的是基础绘图系统的光滑散点图(smoothScatter)。 files_list <- list.files("r-chip/analysis/3-genome-coverage","ReadsCoverage") files_path <- file.path("r-chip/analysis/3-genome-coverage",files_list) input_rep1 <- read.table(files_path[1], sep='\t') input_rep2 <- read.table(files_path[2], sep='\t') input_rep3 <- read.table(files_path[3], sep='\t') chip_rep1 <- read.table(files_path[4], sep='\t') chip_rep2 <- read.table(files_path[5], sep='\t') chip_rep3 <- read.table(files_path[6], sep='\t') pw_plot <- function(x, y, xlab="x", ylab="y", ...){ log2x <- log2(x) log2y <- log2(y) smoothScatter(log2x,log2y, cex=1.2, xlim=c(0,12),ylim=c(0,12), xlab=xlab, ylab=ylab) text(3,10,paste("R = ",round(cor(x,y),2),sep="")) par(mfrow=c(2,3)) pw_plot(chip_rep1[,4], chip_rep2[,4], xlab = "Rep 1 (Log2 Tag Counts)", ylab = "Rep 2 (Log2 Tag Counts)") pw_plot(chip_rep1[,4], chip_rep3[,4], xlab = "Rep 1 (Log2 Tag Counts)", ylab = "Rep 2 (Log2 Tag Counts)") pw_plot(chip_rep2[,4], chip_rep3[,4], xlab = "Rep 1 (Log2 Tag Counts)", ylab = "Rep 2 (Log2 Tag Counts)") pw_plot(input_rep1[,4], input_rep2[,4], xlab = "Rep 1 (Log2 Tag Counts)", ylab = "Rep 2 (Log2 Tag Counts)") pw_plot(input_rep1[,4], input_rep3[,4], xlab = "Rep 1 (Log2 Tag Counts)", ylab = "Rep 3 (Log2 Tag Counts)") pw_plot(input_rep2[,4], input_rep3[,4], xlab = "Rep 2 (Log2 Tag Counts)", ylab = "Rep 3 (Log2 Tag Counts)") BigWig可视化 由于同一个样本间的BAM文件具有很强的相关性,因此可以将这些样本合并起来转换成bigwig格式,这样子在基因组浏览器(例如IGV, UCSC Browser, JBrowse)上方便展示。 如下的代码的目的就是先合并BAM,然后转换成BigWig,拆分成正链和反链进行保存 #!/bin/bash set -e set -u set -o pipefail samples=${1?please provied sample file} threads=${2-8} bs=${3-50} ALIGN_DIR="analysis/2-read-align" BW_DIR="analysis/4-normliazed-bw" mkdir -p $BW_DIR exec 0< $samples cd $ALIGN_DIR while read sample bams=$(ls ${sample}*flt.bam | tr '\n' ' ') if [ ! -f ../$(basename ${BW_DIR})/${sample}.tmp.bam ]; then echo "samtools merge -f -@ ${threads} ../$(basename ${BW_DIR})/${sample}.tmp.bam $bams && samtools index ../$(basename ${BW_DIR})/${sample}.tmp.bam" | bash cd ../../ exec 0< $samples cd ${BW_DIR} while read sample bamCoverage -b ${sample}.tmp.bam -o ${sample}_fwd.bw -of bigwig \ --filterRNAstrand forward --binSize ${bs} --normalizeUsing CPM --effectiveGenomeSize 2864785220 \ --extendReads 150 -p ${threads} 2> ../log/${sample}_fwd.log bamCoverage -b ${sample}.tmp.bam -o ${sample}_rvs.bw -of bigwig \ --filterRNAstrand reverse --binSize ${bs} --normalizeUsing CPM --effectiveGenomeSize 2864785220 \ --extendReads 150 -p ${threads} 2> ../log/${sample}_rvs.log rm -f ${sample}.tmp.bam ${sample}.tmp.bam.bai 得到的BW文件你可以在IGV上初步看看,比如说检查下文章Figure 1(E)提到的CIRH1A基因 #下面将scale等track写入tracklist tracklist<-list() itrack <- IdeogramTrack(genome = "hg19", chromosome = 'chr16',outline=T) tracklist[["itrack"]]<-itrack # 读取BigWig bw_file_path <- "C:/Users/DELL/Desktop/FigureYa/R-ChIP/" bw_file_names <- list.files(bw_file_path, "*.bw") bw_files <- file.path("C:/Users/DELL/Desktop/FigureYa/R-ChIP/", bw_file_names) tracklist[['D210-fwd']] <- DataTrack(range = bw_files[3], genome="hg19", type="histogram", name='D210 + ', ylim=c(0,4), col.histogram="#2167a4", fill.histogram="#2167a4") tracklist[['D210-rvs']] <- DataTrack(range = bw_files[4], genome="hg19", type="histogram", name='D210 - ', ylim=c(4,0), col.histogram="#eb1700", fill.histogram="#eb1700") tracklist[['WKDD-fwd']] <- DataTrack(range = bw_files[9], genome="hg19", type="histogram", name='WKDD + ', ylim=c(0,4), ylab=2, col.histogram="#2167a4", fill.histogram="#2167a4") tracklist[['WKDD-rvs']] <- DataTrack(range = bw_files[10], genome="hg19", type="histogram", name=' WKDD -', ylim=c(4,0), col.histogram="#eb1700", fill.histogram="#eb1700", showAxis=TRUE) #写入比例尺 scalebar <- GenomeAxisTrack(scale=0.25, col="black", fontcolor="black", name="Scale", labelPos="above",showTitle=F) tracklist[["scalebar"]]<-scalebar plotTracks(trackList = tracklist, chromosome = "chr16", from = 69141913, to= 69205033, background.panel = "#f6f6f6", background.title = "white", col.title="black",col.axis="black", rot.title=0,cex.title=0.5,margin=10,title.width=1.75, cex.axis=1 注:这里面有一个希腊字符在不同系统表示有所不同,所以我在复制之后手动修改。 此外,你还需要将里面的(-)和(+)进行替换,因为括号在shell里有特殊含义, 为了保证命名的连贯性,我将(-)替换成-neg,将(+)替换成-pos. sed -i 's/(-)/-neg/; s/(+)/-pos/' sample_name.txt 随后,我们需要从download_table.txt中提取出SRR编号和GSM编号的对应的关系,这个需要用到Linux的文本处理命令 paste <(egrep -o 'GSM[0-9]{6,9}' download_table.txt ) <(egrep -o 'SRR[0-9]{6,9}' download_table.txt) > gsm_srr.txt join <(sort gsm_srr.txt ) <(head -n 24 sample_name.txt | sort) > gsm_srr_sample_name.txt 注: 样本一共有25行,而我们只下载了24个数据,所以删除了最后一行的"K562-WKKD-V5chip" 最后,就是根据生成的文件对样本进行重命名了。 awk '{print "rename " $2 " " $3 " analysis/0-raw-data/" $2 "*.gz"}' gsm_srr_sample_name.txt |bash 这行代码看起来有点复杂,但是其实做的事情就是构建了一系列rename的命令行,然后在bash下执行。 R-ChIP分析 数据预处理 这一部分主要是将序列比对后原来的参考基因组上,标记重复,并且去掉不符合要求的比对。 让我们先写一个比对脚本将序列比对到参考基因组上,脚本命名为03.r_chip_align.sh,存放在scripts下 #!/bin/bash set -e set -u set -o pipefail # configuration threads=8 index=index/hg19 FQ_DIR="analysis/0-raw-data" ALIGN_DIR="analysis/2-read-align" LOG_DIR="analysis/log" TMP_DIR="analysis/tmp" mkdir -p ${ALIGN_DIR} mkdir -p ${LOG_DIR} mkdir -p ${TMP_DIR} samples=${1?missing sample file} exec 0< $samples # alignment while read id; if [ ! -f ${FQ_DIR}/$id.bt2.done ] echo "bowtie2 --very-sensitive-local --mm -p $threads -x $index -U ${FQ_DIR}/$id.fastq.gz 2> ${LOG_DIR}/$id.bt2.log | \ samtools sort -@ 2 -m 1G -T ${TMP_DIR}/${id} -o ${ALIGN_DIR}/${id}.sort.bam \ && touch ${ALIGN_DIR}/$id.bt2.done" | bash 这个脚本的前半部分都在定义各种变量,而#alignment标记的后半部分则是实际运行的比对命令 你会发现我的实际运行部分脚本有点奇怪,我没有直接运行比对,而是用echo通过管道传递给了bash执行。 这样做的原因是, 当我用sed -i 's/| bash/#| bash/ 03.r_chip_align.sh'将| bash替换成#| bash后,运行这个脚本就可以检查代码是否正确,然后再用sed-i s/#| bash/| bash/03.r_chip_align.sh修改回来 。 后面的脚本同理。 接着再写一个脚本用于标记重复04.markdup.sh,也是放在scripts下 #!/bin/bash set -e set -u set -o pipefail # configuration threads=8 index=index/hg19 FQ_DIR="analysis/0-raw-data" ALIGN_DIR="analysis/2-read-align" LOG_DIR="analysis/log" TMP_DIR="analysis/tmp" mkdir -p ${ALIGN_DIR} mkdir -p ${LOG_DIR} mkdir -p ${TMP_DIR} samples=${1?missing sample file} exec 0< $samples # mark duplication while read id; if [ ! -f ${ALIGN_DIR}/${id}.mkdup.done ] echo "sambamba markdup -t $threads ${ALIGN_DIR}/${id}.sort.bam ${ALIGN_DIR}/${id}.mkdup.bam \ && touch ${ALIGN_DIR}/${id}.mkdup.done" | bash 再准备一个脚本用于去除不符合要求的比对,命名为05.bam_filter.sh #!/bin/bash set -e set -u set -o pipefail # configuration threads=8 index=index/hg19 FQ_DIR="analysis/0-raw-data" ALIGN_DIR="analysis/2-read-align" LOG_DIR="analysis/log" TMP_DIR="analysis/tmp" mkdir -p ${ALIGN_DIR} mkdir -p ${LOG_DIR} mkdir -p ${TMP_DIR} samples=${1?missing sample file} exec 0< $samples # filter while read id; if [ ! -f ${ALIGN_DIR}/${id}.flt.done ] echo " samtools view -@ threads -bF 1804 -q 30 ${ALIGN_DIR}/${id}.mkdup.bam -o ${ALIGN_DIR}/${id}.flt.bam\ && samtools index ${ALIGN_DIR}/${id}.flt.bam \ && touch ${ALIGN_DIR}/${id}.flt.done "| bash 最后新建一个samples01.txt用于存放将要处理的样本名 HKE293-D210N-V5ChIP-Rep1 HKE293-D210N-Input-Rep1 HKE293-D210N-V5ChIP-Rep2 HKE293-D210N-Input-Rep2 HKE293-D210N-V5ChIP-Rep3 HKE293-D210N-Input-Rep3 HKE293-WKKD-V5ChIP HKE293-WKKD-Input HKE293-delta-HC-V5ChIP 这样子就可以依次运行脚本了 比对: bash scripts/03.r_chip_align.sh samples01.txt 标记重复:bash scripts/04.bam_markdup.sh samples01.txt 去除不符合要求的比对: bash scripts/05.bam_filter.sh samples01.txt 处理完之后可以对每个样本都进行一次统计,包括如下信息: 处理前的原始reads数 处理后对唯一比对reads数 唯一比对reads数所占原始reads数的比例 这个工作同样可以用shell脚本完成, 脚本为06.sample_align_stat.sh #!/bin/bash set -e set -u set -o pipefail samples=${1?missing sample file} threads=8 ALIGN_DIR="analysis/2-read-align" echo -e "Experiment \t Raw Reads \t Uniquely mapped Reads \t ratio" exec 0< $samples while read sample; total=$( samtools view -@ ${threads} -c ${ALIGN_DIR}/${sample}.sort.bam ) unique=$( samtools view -@ ${threads} -c ${ALIGN_DIR}/${sample}.flt.bam ) ratio=$( echo "scale=2; 100 * $unique / $total " | bc ) echo -e "$sample \t $total \t $unique \t $ratio %" 运行方法是bash scripts/06.sample_align_stat.sh samples01.txt > results/library_stat.txt,运行结果如下,和原本的Table S2对比,你会发现结果基本一致,有出入的地方我推测是标记重复这一步所用软件不同。 Experiment Raw Reads Uniquely mapped Reads ratio HKE293-D210N-V5ChIP-Rep1 22405416 6443979 28.76% HKE293-D210N-Input-Rep1 60302237 25673307 42.57% HKE293-D210N-V5ChIP-Rep2 17763614 11778533 66.30% HKE293-D210N-Input-Rep2 11131443 8553097 76.83% HKE293-D210N-V5ChIP-Rep3 8799855 5640375 64.09% HKE293-D210N-Input-Rep3 4529910 3209275 70.84% HKE293-WKKD-V5ChIP 12734577 8612940 67.63% HKE293-WKKD-Input 8830478 6643507 75.23% HKE293-delta-HC-V5ChIP 25174573 9252009 36.75%

R-loop数据分析之R-ChIP(环境准备)

提高自己分析能力的一个好的方法就是重复别人文章里的分析策略,所以这里会尝试对第一篇介绍R-ChIP技术文章"R-ChIP Using Inactive RNase H Reveals Dynamic Coupling of R-loops with Transcriptional Pausing at Gene Promoters"里的所有分析进行重复,我重复所用代码会更新在我的GitHub上,地址为https://github.com/xuzhougeng/R-ChIP-data-analysis 选择这篇文章进行重复的理由有三点: 一:最近要探索R-loop数据分析流程 二:这篇文章的通讯作者是大牛,Xiang-Dong Fu 三:这篇文章将分析所用代码都托管在https://github.com/Jia-Yu-Chen 我整理下和数据分析有关的几个知识点: R-loop是一种RNA/DNA三链结构体,与基因组稳定性和转录调控有关。 通过电镜观察,R-loop大小在150~500bp之间。 硫酸氢盐测序(bisulfate sequencing)表明R-loop主要出现在基因启动子的下游。 R-loop所在非模板链(又称编码链)具有很强的序列偏好性,计算方式为(G-C)/(G+C) R-loop的高通量分析方法目前都是依赖于S9.6抗体捕获RNA/DNA杂合体,然后超声打断或酶切,如果后续对DNA进行测序,那就是DRIP-seq(DNA:RNA immunoprecipitation [DRIP] sequencing),如果后续对RNA逆转成的cDNA继续测序,那就是 [DRIPc]-seq(DNA:RNA immunoprecipitation followed by cDNA conversion)。 然而酶切的分辨率不够,超声又容易破坏脆弱的R-loop结构,于是就导致目前很多文献报道有矛盾。 这篇文章就开发了一种新方法,基于RNase H的体内R-loop谱检测策略。作者构建一种没有催化活性,且在C端有一个V5标签的RNASE H1,RNASEH1与RNA/DNA结合,超声打碎,用anti-V5抗体进行染色体免疫共沉淀(ChIP)。随后RNA/DNA杂合体转换成双链DNA(ds-DNA), 之后便是链特异性测序。 关于链特异性测序,推荐拜读链特异性测序那点事 文章中"Software and Algorithms"这部分列出了分析主要所用的软件,加上下载SRA数据所需工具和一些常用软件,一共要安装的软件如下: SRA Toolkit: 数据下载工具 Bowtie2: 比对工具 SAMtools: SAM格式处理工具 BEDtools: BED格式处理工具 MACS2: 比对后找peak R: 统计作图 Ngsplot: 可视化工具 Deeptools: BAM文件分析工具, 可作图。 软件安装部分此处不介绍,毕竟如果你连软件安装都有困难,那你应该需要先学点Linux基础,或者去看生信必修课之软件安装 分析项目搭建 使用mkdir创建项目文件夹,用于存放后续分析的所用到的数据、中间文件和结果 mkdir -p r-chip/{analysis/0-raw-data,index,scripts,results} 个人习惯,在项目根目录下创建了四个文件夹 analysis: 存放原始数据、中间文件 index: 存放比对软件索引 scripts: 存放分析中用到的脚本 results: 存放可用于放在文章中的结果 后续所有的操作都默认在r-chip下进行,除非特别说明。 根据文章提供的GEO编号(GEO: GSE97072)在NCBI上检索, 按照如下步骤获取该编号下所有数据的元信息, 我将其重命名为"download_table.txt"然后上传到服务器, 。 使用如下命令进行数据下载 tail -n+2 download_table.txt | cut -f 6 | xargs -i prefetch {} >> download.log & 下载的数据默认情况下存放在~/ncbi/public/sra, 需要用fastq-dump解压缩到analysis/0-raw-data. fastq-dump的使用说明见Fastq-dump: 一个神奇的软件 新建一个脚本,叫做uncompress.sh,存放在scripts文件下,代码如下 #!/bin/bash set -e set -o pipefail set -u tail -n+2 download_table.txt | cut -f 6 | while read id; fastq-dump --gzip --split-3 --defline-qual '+' --defline-seq '@$ac-$si/$ri' &id -O analysis/0-raw-data & 然后用bash scripts/uncompress.sh运行。 注意:这是单端测序,所以每个SRR只会解压缩出一个文件 此外还需要下载human genome (hg19)的bowtie2索引,用于后续bowtie2比对。 curl -s ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/hg19.zip -o index/hg19.zip & cd index unzip hg19.zip

作图的目的是希望在图里面发现问题或者解释问题,当然更本质一点就是你想解决什么问题? 前几天做了一个PCA的图,图是画出来了,但是问题有很多,比如说主成分是是啥意思,图里面的箭头有什么含义?为了不做无意义的重复,所以写一篇文章尝试做一个解释。 我们以R语言自带的数据集iris作为例子来演示。 data(iris) iris翻译成中文就是鸢(yuan, 第一声)尾花(如下图), 我建议你在R语言里用?iris了解更多这个数据集的出处。 先让我们用head简单看下这个数据集的前面几行。其中前面四列是一朵鸢尾花的一些形态特征衡量值,自行翻译各个单词的含义。最后一列是这朵花的所属物种,根据分类,这些花来自于"setosa, versicolor,virginica" > head(iris) Sepal.Length Sepal.Width Petal.Length Petal.Width Species 1 5.1 3.5 1.4 0.2 setosa 2 4.9 3.0 1.4 0.2 setosa 3 4.7 3.2 1.3 0.2 setosa 4 4.6 3.1 1.5 0.2 setosa 5 5.0 3.6 1.4 0.2 setosa 6 5.4 3.9 1.7 0.4 setosa 在平常人眼里, 看到这些数值可能没有任何想法,也不知道这到底是什么品种,但是对于一个研究鸢尾花的人,这些数值可能会在他们的脑海里重构出一朵鸢尾花,然后迅速判断出这是什么品种。我和阅读文章的人都一样,也不知道如何区分鸢尾花的品种,如果想去研究这种规律,应该怎么做呢? 人类的学习过程过程大多数多看多提出一些模型规律,然后基于这个规律去判断,如果错了改正模型,直到找到一个好用的标准为止。因此,后面我们也就是在当前的数据集下多试试,看看能不能提出一些模型。 可视化是非常强大的工具,让我们先看看能不能用“ Sepal.Length”和“ Sepal.Width”进行区分 library(ggplot2) ggplot(iris,aes(x=Sepal.Length, y=Sepal.Width)) + geom_point(aes(colour=Species)) 从上图中,我们发现根据萼片(sepal)的宽度和长度似乎能够区分出"setosa"和"versicolor","virginica",但是"versicolor","virginica"不太好分。 让我们再试试花瓣(petal)的长度和宽度。从下图,我们可以发现在这两个属性下,不同品种的鸢尾花似乎有了明显的分界线。 ggplot(iris,aes(x=Petal.Length, y=Petal.Width)) + geom_point(aes(colour=Species)) 当花瓣长度小于2时就一定是setosa 当花瓣的宽度在0.75~1.75之间,花瓣的长度在3~5,则是versicolor, 当花瓣的宽度大于1.75时,花瓣的长度大于5,则是virginica 这个模型依旧还不完善,存在一些点难以区分。不过我们还可以以花瓣长度和萼片长度,花瓣宽度和萼片宽度等你能想到的组合进行作图,说不定能找到更好的模型。更好的情况是,如果人类能够想象出一个四维的空间,在这个空间里面同时展现这四个变量,或许就能更加直观的找到规律,可惜,我们目前还没有这种能力 有没有一种方法能够把这种多维数据降维,也就是抓住数据集的主要矛盾呢?当然是有,这个方法就是100多年前,由卡尔·皮尔森提出的 主成分分析,通过将原来众多具有一定相关性的变量,重新组合成一组新的互相无关的综合变量。比如说我们可以将花瓣长度和宽度组合成花瓣的面积,当然实际计算可能更加复杂,你需要看机器学习中的数学(4)-线性判别分析(LDA), 主成分分析(PCA)。 反正PCA最终的目的就是让原本的数据在降维后方差最大,也就是能把数据分开。我们一般人记住这一点,然后会用软件,能解读结果就好。 R语言能做主成分分析功能很多,比如说prcomp,princomp,principal. 这里用到是prcomp,它使用的是奇异值分解而不是特征值分解,至于这两个有啥区别,其实这不是我关心的,当然你想深入了解到话,建议阅读机器学习中的数学(5)-强大的矩阵奇异值分解(SVD)及其应用。 代码如下: pca.results <- prcomp(iris[1:4],center = TRUE, retx = T) prcomp会返回一个列表,里面包含如下内容,你可以用print(pca.results)查看 sdev: 对应每个主成分的标准偏差,值越大说明解释变异越大。 rotation: 变量的负荷(loading)矩阵,用于衡量一个变量在各个主成分重要性。 x: 原矩阵经线性变化后的矩阵(下图不显示,可以用pca.results$x提取) 虽然有现成的工具提取pca.results里的数据进行作图,但这样无助于理解,我们需要自己手动提取结果进行作图。 iris_rotate_df <- as.data.frame(pca.results$x[,c("PC1","PC2")]) iris_rotate_df$Species <- iris$Species p <- ggplot(iris_rotate_df, aes(x=PC1, y=PC2)) + geom_point(aes(colour=Species)) print(p) 这个结果虽然也无法完美地区分出"versicolor"和"virginica",但是比我们盲目的寻找好多了。这里的主成分是原来的各个变量线性组合结果。 当然我们还有一个问题,各个变量在不同的主成分里的权重是多少?这个就要看负荷矩阵了,也就是Rotation部分。 在PCA中,我们将协方差矩阵分成了标量部分(特征值)和有向部分(特征向量),通过计算特征向量和开方后特征值的乘积,就称之为负荷(loading) 我们发现在PC1里面Petal.Length的绝对值最大,而在PC2里面,Sepal.Width的绝对值最大,于是我们就可以从原来的数据集中提取这两个变量进行作图了。 ggplot(iris,aes(x=Sepal.Width, y=Petal.Length)) + geom_point(aes(colour=Species)) 下面的作图直接画图的结果,是不是感觉和上面的PCA图很像。当我手动将其旋转后得到下面的右图,你会发现这和PCA的结果几乎一摸一样。 那么我们如何在图上展示各个变量在各个PC的占比呢?这个就需要画几个箭头了, p + annotate("segment", x=0,xend=0.35,y=0,yend=-0.65,arrow=arrow()) + annotate("text", x=0.35,y=-0.68, label="Sepal.Length") + annotate("segment", x=0,xend=-0.08,y=0,yend=-0.73,arrow=arrow()) + annotate("text", x=-0.08,y=-0.75, label="Sepal.Width") 下图中"Septal.Width"朝下,Petal.Length朝右,如果你手动旋转一下,就是差不多是以"Septal.With"为X轴,"Petal.Length"为Y轴的结果了。 所以箭头的朝向不重要,重点是长度。

我终于知道为啥"find /home -name .bashrc > list 2>list"是错误的啦

在《鸟哥的Linux私房菜》第三版的329-330,鸟哥介绍了数据流的重定向,里面提到一点,如何你想把正确与错误数据都写到一个文件中,这个时候用下面第一行代码是错误的,而使用第二行和第三行代码才是正确的 find /home -name .bashrc > list 2>list find /home -name .bashrc > list 2>&1 find /home -name .bashrc &> list 至于第一行错误的原因,鸟哥说这是因为两条数据同时写入一个文件,有没有使用特殊的语法,此时两条数据可能会交叉写入该文件内,造成次序的错乱。 在刚开始学Linux的时候,我只记住了结论,就是用&> 将标准输出和标准错误输出写入到一个文件中,而2>&1要记忆4个字符,还得记住要放在文件的后面,脑容量不够,就没有记住。而这两天读到了《Linux命令行与shell脚本编程大全》第3版第15章的时候,我终于明白了为什么第一行代码会出现次序错乱,2>&1为什么这些写是正确的。 首先,你需要理解1>和2>里的数字1和2指的是文件描述符(file descriptor)。 所谓的文件描述符是一个非负整数,用于标识打开的文件。在bash shell中,1表示STDOUT(标准输出),2表示STDERR(标准错误), 当然还有一个0表示的STDIN(标准输入)。我们还可以定义其他整数,但这个就比较复杂,算是高级技巧了。 其次,你要知道">"和另一个符号"<"是两个常用重定向符号。所谓的重定向,就是改变数据流的原本走向,本来输出到屏幕的输出写入到文件中,本来需要从文件读取的数据,可以从键盘上输入 那么&符号是什么含义呢?在C语言里&是一个求址符号,因此&1你可以理解成获取标准输出对象的地址,而>&1就表示把输入定向到了标准输出, >&2也就是把输出定向到了标准错误输出。那么我们就可以在脚本里面写一些标准错误输出了, 比如创建一个脚本, test.sh #!/bin/bash echo "This is a first normal output" echo "This is a first error" >&2 echo "This is a second error" >&2 echo "This is a second nromal output" 以bash test.sh 1> a 2>a运行的话,由于同时对a文件进行写入操作,不同的磁盘写入速度,会有不同的结果,因为前面的输出很可能会后面的覆盖了。 cat a This is a first error This is This is a second nromal output 从我的输出结果可以推测出,首先是输出“This is a first normal output”,但是被后来者“This is a first error”给覆盖了。然后“This is a second error"写到一半,就被“This is a second nromal output”给覆盖了。 而以bash test.sh 1> a 2>&1的话,就会把标准错误输出定向到标准输入上,就会以此输出了。 This is a first normal output This is a first error This is a second error This is a second nromal output 考虑到这个操作比较常用,每次都要写2>&1太麻烦了,于是就定义了一个&>符号方便使用。

使用MAKER进行基因注释(高级篇之SNAP模型训练)

训练 ab initio 基因预测工具(以SNAP为例) 对于一个新的物种而言,你大概率是没有一个高质量的基因模型去进行基因预测。但是我们可以利用EST序列(少部分物种估计有)、二代测序数据、同源物种蛋白序列,先直接用Maker做基因注释,尽管得到的模型可能不是特别的完美,但可以作为输入反复迭代运行Maker,从而提高最终的表现。 这次使用的是下载的练习数据集(见附录) cd ~/maker_tutorial/example_02_abinitio 同样,让我们先构建配置文件,并修改如下配置 maker -CTL vim maker_opts.ctl # modify the following line genome=pyu_contig.fasta est=pyu_est.fasta protein=sp_protein.fasta est2genome=1 protein2genome=1 这里的"est2genome"和"protein2genome"表示直接从EST序列和同源但序列中推测基因结构,当然这肯定不靠谱。不过没有关系,我们的目标是将其作为输入用于训练而已。 运行预测程序,大约需要20分钟 ~/opt/biosoft/maker/bin/maker &> maker.log & 那么下一步就是收集所有的GFF文件,整理成SNAP所需的ZFF格式 mkdir snap cd snap ~/opt/biosoft/maker/bin/gff3_merge -d ../pyu_contig1.maker.output/pyu_contig1_master_datastore_index.log ~/opt/biosoft/maker/bin/maker2zff pyu_contig1.all.gff 于是我们就会在snap文件下得到"genome.ann"和"genome.dna". 在这两个文件的基础上,我们就可以参考SNAP的文档开始训练 可以先用fathom genome.ann genome.dna -gene-stats了解基因的一些信息,比如说这里的测试数据集就有153个基因,几乎平均的分布在正负链上。 1 sequences 0.525725 avg GC fraction (min=0.525725 max=0.525725) 153 genes (plus=79 minus=74) 5 (0.032680) single-exon 148 (0.967320) multi-exon 130.782104 mean exon (min=3 max=704) 87.851593 mean intron (min=61 max=384) 此外还可以用fathom genome.ann genome.dna -validate检查下是否有明显的错误,这里的153个基因有106个warning,警告类型粗略看了一眼基本都是CDS不完整。 后续就可以开始参数预测。步骤是,先用fathom -genome.ann genome.dna -categorize 1000将序列分类,这里的1000表示基因两侧会额外有1000bp的序列。该参数推荐使用基因一半的长度,如果基因比较稠密则要调低。这一步会生成如下文件: alt.ann, alt.dna (genes with alternative splicing) err.ann, err.dna (genes that have errors) olp.ann, olp.dna (genes that overlap other genes) wrn.ann, wrn.dna (genes with warnings) uni.ann, uni.dna (single gene per sequence) 这里只用第四类的基因,也就是每个序列上只有一个基因。用fathom uni.ann uni.dna -export 1000 -plus只输出unigene中正链基因,这一步同样会生成四个文件 export.aa 每个基因的蛋白序列 export.ann 正链的基因结构 export.dna 正链的DNA export.tx 每个基因的转录本 接着让forge负责预测参数, 由于输出会很多,所以建议创建个文件夹 mkdir params cd params forge ../export.ann ../export.dna cd .. 最后是hmm-assembler.pl构建HMM,即基因模型文件, hmm-assembler pyu params > pyu1.hmm 完成SNAP的模型构建后,修改"maker_opts.ctl"用以增加该文件,并不再用est和protein直接推测基因结构。 snaphmm=pyu1.hmm est2genome=0 protein2genome=0 再一次运行maker ~/opt/biosoft/maker/bin/maker &> maker.log & 这次结果会比上一次有很明显的提升,你可以重复上面的代码从而进一步提高SNAP的模型。

集群任务管理系统SGE的简明教程

我用的一个服务器上装了一个集群管理工具(SGE, Sun Grid Engine), 用于从登陆节点上向计算节点进行任务投递。一开始,不太会用,但是经过一段时间的摸索学习后,终于能顺手的用起来了。 在使用SGE之前,你得先了解SGE到底做了什么事情. SGE或者其他集群管理工作做的事情就是将用户投递的任务进行排队,然后将任务交给能够运行的结算节点执行,工作流程可以分为四步: 接受用户投放的任务 在任务运行以前,将任务放到一个存储区域 发送任务到一个执行设备,并监控任务的运行 运行结束写回结果并记录运行日志 SGE的常用命令 SGE中投递任务所用到的命令是qsub. 最简单的用法是下面这种,即,将要执行的命令通过标准输入的方式传递给qsub echo "ls -l " | qsub 投递之后可以用qstat查看自己投递的任务的运行情况,如下图 任务投递情况 第一列是任务编号, 第二列是优先级,第三列是任务名字,在参数里没有特别说明的情况下,SGE会用任务的来源进行命令,STDIN表示来自于标准输入,第四列是用户名,第五列是运行状态("r"表示运行中), 第六列表示任务投递和开始时间,第七列是任务投递的节点,第8列则是要申请的线程数。在执行完成后会在家目录下生成"STDIN.e7883"和"STDIN.o7883", 其中7883就是任务编号, 前者存放标准错误输出, 后者存放标准输出, 因此"cat STDIN.o7883"的内容就是ls -l的内容。 另一种方法是先写一个脚本然后投递,比如先编辑一个文件"ls.sh", 内容如下,然后用"qsub ls.sh"投递任务。 ls -l 跟之前一样,最后在家目录下产生了"ls.sh.exxxx"和"ls.sh.exxxx"两个文件 当然实际时肯定没有那么简单,我们需要增加各种参数来调整qsub的行为,用qsub -help可以看完整的参数,但是常用的为如下几个 -q xxx : 指定要投递到的队列,如果不指定的话,SGE会在用户可使用的队列中选择一个满足要求的队列。 -V: 将当前的环境变量传递到执行命令的节点中 cwd: 在当前目录下执行任务, sge的日志会输出到当前路径。 不增加该指令,所有投递的任务都会在家目录下执行 -l resource=value: 请求资源数, 例如 -l vf=25G -l h=node1 就是任务的预估内存要25G(内存估计的值应稍微大于真实的内存,内存预估偏小可能会导致节点跑挂), 申请在node1上运行 -S /bin/bash: 表示在bash环境下执行命令。默认tcsh. -pe openmpi 4: 表示使用openmpi进行并行运算,且申请的线程是4, -N 任务名: 手动执行任务的名字 -j y|n :是否将标准输入和标准输入合并成一个文件 -sync y|n: 是否等待任务结束,返回退出码 -o path: 指定标准输出的文件夹 那么接下来就可以添加这些参数运行一些命令了,例如在命令行里投递一个比对任务 echo "bowtie2 -p 8 -x index/ref -1 data/A_1.fq -2 data/A_2.fq | samtools sort > A.bam" | qsub -V -cwd -l vf=25G -S /bin/bash -pe openmpi 8 -N A.bt2 这些参数除了在外部设置外,还可以在shell脚本里设置,如下 #!/bin/bash #$ -S /bin/bash #$ -V #$ -cwd #$ -l vf=25G #$ -pe openmpi 8 #$ -N a.bt2 bowtie2 -p 8 -x index/ref -1 data/A_1.fq -2 data/A_2.fq | samtools sort > A.bam 除了任务投递外,查询任务也是一个非常常用的命令,除了刚才直接用qstat查看,还有如下参数比较好用 qstat -f # 查看用户任务 qstat -j jobId # 按任务id查看 qstat -explain a|c|A|E -j jobID # 查看任务任务并给出解释 qstat -u user # 按用户查看 任务状态: qw: 表示等待状态 hqw: 任务挂起等待中,待依赖的任务完成后执行 Eqw: 投递任务出错 r: 表示任务正在运行 s: 暂时挂起 dr: 节点挂了之后,删除任务就会出现这个状态,只有节点重启之后,任务才会消失 任务删除也比较重要,毕竟偶尔会出现任务投递出错的情况 qdel -j 1111 删除任务号为1111的任务 qrsh:与qsub相比,是交互式的投递任务,注意参数:-now yes|no默认设置为yes 若设置为yes,立即调度作业,如果没有可用资源,则拒绝作业,任务投递失败,任务状态为Eqw。 若设置为no,调度时如果没有可用资源,则将作业排入队列,等待调度。 例子: qrsh -l vf=*G -q all.q -now no -w n *sh qacct 从集群日志中抽取任意账户信息 qalter 更改已提交但正处于暂挂状态的作业的属性 qconf 为集群和队列配置提供用户界面 qconf -spl查看可用并行环境 qhold 阻止已提交作业的执行 qhost 显示SGE执行主机(即各个计算节点)的状态信息 qhost -j按照节点显示任务 qhost -F展示每个节点的资源 qlogin 启动telnet或类似的登录会话。 案例:一个投递比对任务的简单脚本 #!/bin/bash set -e set -u set -o pipefail threads=8 index=index/hg19 FQ_DIR="analysis/0-raw-data" ALIGN_DIR="analysis/2-read-align" LOG_DIR="analysis/log" TMP_DIR="analysis/tmp" mkdir -p ${ALIGN_DIR} mkdir -p ${LOG_DIR} mkdir -p ${TMP_DIR} tail -n+2 download_table.txt | cut -f 6 | while read id; echo " bowtie2 --very-sensitive-local --mm -p $threads -x $index -U ${FQ_DIR}/$id.fastq.gz 2> ${LOG_DIR}/$id.bt2.log | \ samtools sort -@ 2 -m 1G -T ${TMP_DIR}/${id} -o ${ALIGN_DIR}/${id}.sort.bam" | qsub -V -cwd -pe openmpi $threads -N ${id}_bt2 -q all.q -S /bin/bash

基因组组装完成后,或者是完成了草图,就不可避免遇到一个问题,需要对基因组序列进行注释。注释之前首先得构建基因模型,有三种策略: 从头注释(de novo prediction):通过已有的概率模型来预测基因结构,在预测剪切位点和UTR区准确性较低 同源预测(homology-based prediction):有一些基因蛋白在相近物种间的保守型搞,所以可以使用已有的高质量近缘物种注释信息通过序列联配的方式确定外显子边界和剪切位点 基于转录组预测(transcriptome-based prediction):通过物种的RNA-seq数据辅助注释,能够较为准确的确定剪切位点和外显子区域。 每一种方法都有自己的优缺点,所以最后需要用EvidenceModeler(EVM)和GLEAN工具进行整合,合并成完整的基因结构。基于可靠的基因结构,后续可才是功能注释,蛋白功能域注释,基因本体论注释,通路注释等。 那么基因注释重要吗?可以说是非常重要了,尤其是高通量测序非常便宜的现在。你可以花不到一万的价格对600M的物种进行100X的普通文库测序,然后拼接出草图。但是这个草图的价值还需要你进行注释后才能显现出来。有可能你和诺贝尔奖就差一个注释的基因组。 从案例中学习套路 陆地棉基因组注释 文章标题为“Sequencing of allotetraploid cotton (Gossypium hirsutum L. acc. TM-1) provides a resource for fiber improvement”. 同源注释:从Phytozome上下载了7个植物的基因组蛋白序列(Arabidopsis thaliana, Carica papaya, Glycine max, G. raimondii, Populus trichocarpa, Theobroma cacao and Vitis vinifera), 使用 TblastN 将蛋白序列比对到组装序列上,E-value的阈值为1e-5. 将不同蛋白的BLAST的hits用 Solar 软件进行合并。GeneWise 根据每个BLAST hit的对应基因区域预测完整的基因结构。 从头预测:先得构建repeat-mask genome, 在这个基础上就用 August, Genescan, GlimmerHMM, Geneid 和 SNAP 预测编码区 转录组预测:用Tophat将RNA-seq数据比对到组装序列上,然后用cufflinks组装转录本形成基因模型。 综上,使用 EvidenceModeler(EVM) 将上面的结果组装成非冗余的基因结构。进一步根据Cscore > 0.5,peptide coverage > 0.5 和CDS overlaping with TE进行筛选。还有过滤掉超过30%编码区被Pfam或Interprot TE domain的注释的基因模型。 这些基因模型使用BLASTP进行功能注释,所用数据库为SWiss-Prot和TrEMBL.蛋白功能使用InterProScan和HMMER注释,数据库为InterPro和Pfam。GO注释则是直接雇佣InterPro和Pfam注释得到的对应entry。通路注释使用KEGG数据库。 Cardamine hirsuta基因组注释 文章标题为“The Cardamine hirsuta genome offers insight into the evolution of morphological diversity”。 同源注释:使用 GenomeThreader 以拟南芥为剪切模型,以及PlantsGDB resourc上 Brassica rapa (v1.1), A. thaliana(TAIR10), A. lyrata (v6), tomato (v3.6), poplar (v2) 和 A. thaliana (version PUT-169), B. napus (version PUT-172) EST assemblies 的完整的代表性蛋白集。 转录本预测: 将 C. hirsuta RNA-seq数据比对到基因序列,然后用cufflinks拼接 从头预测:转录本预测得到的潜在蛋白编码转录本使用网页工具 ORFpredictor 进行预测, 同时用 blastx 和 A. thalina 进行比较,选择90%序列相似度和最高5%长度差异的部分从而保证保留完整的编码框(有启动子和终止子)。 这些基因模型根据相互之间的相似度和重叠度进行聚类,高度相似(>95)从聚类中剔除,保证非冗余训练集。为了训练gene finder, 它们选随机选取了2000个位点,20%是单个外显子基因。从头预测工具为 August , GlimmerHMM, Geneid 和 SNAP . 此外还用了Fgenesh+, 以双子叶特异矩阵为参数进行预测。 最后使用JIGSAW算法根据以上结果进行训练,随后再次用JIGSAW对每个基因模型计算统计学权重。 可变剪切模型则是基于苗、叶、花和果实的RNA-seq比对组装结果。 GO注释使用AHRD流程 举的2个例子都是植物,主要是植物基因组不仅是组装,注释都是一大难题。因为植物基因组有大量的重复区,假基因,还有很多新的蛋白编码基因和非编码基因,比如说玉米基因组80%以上都是重复区域。然后当我检索这两篇文章所用工具的时候,我不经意或者说不可避免就遇到了这个网站 http://www.plantgdb.org/ , 一个整合植物基因组学工具和资源的网站,但是这个网站似乎2年没有更新了。当然这个网站也挺不错,http://bioservices.usd.edu/gsap.html, 他给出了一套完整的注释流程以及每一步的输入和输出情况。 此外,2017年在《Briefings in Bioinformatics》发表的"Plant genome and transcriptome annotations: from misconceptions to simple solution" 则是从五个角度对植物基因组注释做了很完整的总结 植物科学的常见本体 功能注释的常用数据库和资源 已注释的植物基因组意味着什么 一个自动化注释流程 一个参考流程图,用来说明使用公用数据库注释植物基因组/转录组的常规步骤 gene structure 在正式启动基因组注释项目之前,需要先检查组装是否合格,比如contig N50的长度是否大于基因的平均长度,使用BUSCO/CEGMA检查基因的完整性,如果不满足要求,可能输出结果中大部分的contig中都不存在一个完整的基因结构。当组装得到的contig符合要求时,就可以开始基因组注释环节,这一步分为三步:基因结构预测,基因功能注释,可视化和质控。 基因组结构注释 基因结构注释应是功能注释的先决条件,完整的真核生物基因组注释流程需要如下步骤: 必要的基因组重复序列屏蔽 从头寻找基因, 可用工具为: GeneMarkHMM, FGENESH, Augustus, SNAP, GlimmerHMM, Genscan 同源蛋白预测, 内含子分析: GeneWIse, Exonerate, GenomeThreader 将EST序列,全长cDNA序列和Trinity/Cufflinks/Stringtie组装的转录组和基因组联配 如果第4步用到了多个数据来源,使用PASA基于重叠情况进行联配 使用EvidenceModler根据上述结果进行整合 使用PASA更新EVM的一致性预测,增加UTR注释和可变剪切注释 必要的人工检查 基本上是套路化的分析流程,也就有一些工具通过整合几步开发了流程管理工具,比如说BRAKER结合了GeneMark和Augustus,MAKER2整合了SNAP,Exonerate,虽然BRAKER说自己的效果比MAKER2好,但是用的人似乎不多,根据web of knowledge统计,两者的引用率分别是44,283, 当然BRAKER是2016,MAKER2是2011,后者在时间上有优势。 这里准备先按部就班的按照流程进行注释,所用的数据是 Cardamine hirsuta , 数据下载方式如下 # Cardamine hirsutat基因组数据 mkdir chi_annotation && cd chi_annotation wget http://chi.mpipz.mpg.de/download/sequences/chi_v1.fa cat chi_v1.fa | tr 'atcg' 'ATCG' > chi_unmasked.fa # 注释结果 wget http://chi.mpipz.mpg.de/download/annotations/carhr38.gff # Cardamine hirsutat转录组数据 mkdir rna-seq && cd rna-seq wget -4 -q -A '*.fastq.gz' -np -nd -r 2 http://chi.mpipz.mpg.de/download/fruit_rnaseq/cardamine_hirsuta/ & wget -4 -q -A '*.fastq.gz' -np -nd -r 2 http://chi.mpipz.mpg.de/download/leaf_rnaseq/cardamine_hirsuta/ & 软件安装不在正文中出现,会放在附录中,除了某些特别复杂的软件。 01-重复序列屏蔽 重复屏蔽:真核生物的基因组存在大量的重复序列,植物基因组的重复序列甚至可以高达80%。尽管重复序列对维持染色体的空间结构、基因的表达调控、遗传重组等都具有重要作用,但是却会导致BLAST的结果出现大量假阳性,增加基因结构的预测的计算压力甚至影响注释正确性。基因组中的重复按照序列特征可以分为两类:串联重复(tandem repeats)和散在重复(interspersed repeats). 人类中的重复序列划分 鉴定基因组重复区域的方法有两种:一种基于文库(library)的同源(homology)方法,该文库收集了其他物种的某一种重复的一致性序列,通过相似性来鉴定重复;另一种是从头预测(de novo),将序列和自己比较或者是高频K-mer来鉴定重复。 目前重复序列注释主要软件就是RepeatMasker和RepeatModel。这里要注意分析的fasta的ID不能过长,不然会报错。如果序列ID过长可以使用bioawk进行转换,后续用到RepatModel不支持多行存放序列的fasta格式。 直接使用同源注释工具RepeatMasker寻找重复序列: mkdir 00-RepeatMask ~/opt/biosoft/RepeatMasker/RepeatMasker -e ncbi -species arabidopsis -pa 40 -gff -dir 00-RepeatMask/ chi_unmasked.fa # -e ncbi # -species 选择物种 用~/opt/biosoft/RepeatMasker/util/queryRepeatDatabase.pl -tree 了解 # -lib 增加额外数据库, # -pa 并行计算 # -gff 输出gff注释 # -dir 输出路径 # annotation with the library produced by RepeatModel 输出结果中主要关注如下三个(其中xxx表示一类文件名) xxx.fa.masked, 将重复序列用N代替 xxx.fa.out.gff, 以gff2形式存放重复序列出现的位置 xxx.fa.tbl, 该文件记录着分类信息 cat 00-RepeatMask/chi_unmasked.fa.tbl ================================================== file name: chi_unmasked.fa sequences: 624 total length: 198654690 bp (191241357 bp excl N/X-runs) GC level: 35.24 % bases masked: 35410625 bp ( 17.83 %) ================================================== 也就是说该物种198M中有将近18%的重复序列,作为参考,拟南芥125Mb 14%重复序列, 水稻389M,36%重复,人类基因组是3G,50%左右的重复序列。 使用最后的chi_unmasked.fa.masked用于下一步的基因结构预测。 注:当然也可以用RepeatModel进行从头预测,得到的预测结果后续可以整合到RepeatMasker # de novo predict ~/opt/biosoft/RepeatModeler-open-1.0.11/BuildDatabase -name test -engine ncbi output.fa ~/opt/biosoft/RepeatModeler-open-1.0.11/RepeatModeler -database test 这一步速度极其慢,由于我们的目的只是获取屏蔽后序列降低后续从头预测的压力,所以可以先不做这一步。在后续分析重复序列在基因组进化上的作用时可以做这一步。下 如果从头预测的结果与同源预测的结果有30%以上的overlap,并且分类不一致,会把从头预测的结果过滤掉。从头预测与同源预测结果有overlap,但是分类一致的,都会保留。但是统计的时候不会重复统计。 02-从头(ab initio)预测基因 基于已有模型或无监督训练 目前的从头预测软件大多是基于HMM(隐马尔科夫链)和贝叶斯理论,通过已有物种的注释信息对软件进行训练,从训练结果中去推断一段基因序列中可能的结构,在这方面做的最好的工具是AUGUSTUS 它可以仅使用序列信息进行预测,也可以整合EST, cDNA, RNA-seq数据作为先验模型进行预测。 AUGUSTUS的无root安装比较麻烦,我折腾了好几天最后卒,不过辛亏有bioconda,conda create -n annotation augustus=3.3. 它的使用看起来很简单,我们可以尝试使用一段拟南芥已知的基因序列让其预测,比如前8k序列 seqkit faidx TAIR10.fa Chr1:1-8000 > test.fa augustus --speices=arabidopsis test.fa > test.gff 如果仅仅看两者的CDS区,结果完全一致,相当于看过一遍参考答案去做题目,题目都做对了。 注:已经被训练的物种信息可以用augustus --species=help查看。 在不使用RNA-seq数据的情况下,可以基于拟南芥的训练模型进行预测,采用下面的方式多条染色体并行augustus mkdir 01-augustsus && cd 01-augustsus ln ../00-RepeatMask/chi_unmasked.fa.masked genome.fa seqkit split genome.fa #结果文件在genome.fa.split find genome.fa.split/ -type f -name "*.fa" | parallel -j 30 augustus --species=arabidopsis --gff3=on >> temp.gff #并行处理 join_aug_pred.pl < temp.gff | grep -v '^#' > temp.joined.gff bedtools sort -i temp.joined.gff > augustsus.gff AUGUSTUS依赖于已有的模型,而GeneMark-ES/ET则是唯一一款支持无监督训练模型,之后再识别真核基因组蛋白编码区的工具。 gmes_petap.pl --ES --sequence genome.fa --cores 50 最后得到的是genemark.gtf,是标准的GTF格式,可以使用Sequence Ontology Project提供的gtf2gff3.pl进行转换 wget http://genes.mit.edu/burgelab/miso/scripts/gtf2gff3.pl chmod 755 gtf2gff3.pl gtf2gff3.pl genemark.gtf | bedtools sort -i - > genemark.gff 不同从头预测软件的实际效果可以通过在IGV中加载文章提供的gff文件和预测后的gff文件进行比较,一般会存在如下几个问题: 基因多了,或者少了,也就是假阳性和假阴性现象 UTR区域难以预测,这个比较正常 未正确识别可变剪切位点,导致前后几个基因识别成一个基因 考虑到转录组测序已经非常便宜,可以通过该物种的RNA-seq提供覆盖度信息进行预测。 基于转录组数据预测 根据已有的模型或者自训练可以正确预测很大一部分的基因,但如果需要提高预测的正确性,还需要额外的信息。在过去就需要提供物种本身的cDNA, EST,而现在更多的是基于转录组序列进行训练。尽管RNA-seq数据在基因组上的比对情况能够推测出内含子位置,根据覆盖度可以推测出外显子和非编码区的边界,但是仅仅依赖于RNA-seq的覆盖不能可信地推测出蛋白编码区(Hoff K.J. Stanke M. 2015). AUGUSTUS可以利用转录组比对数据中的位置信息来训练模型,GeneMark-ET可以利用RNA-seq得到的内含子位点信息自我训练HMM参数,进行基因预测。BRAKER2将两者进行整合,使用GeneMark-ET根据RNA-seq无监督训练模型寻找基因,然后用AUGUSTUS进行模型训练,最后完成基因预测 mkdir index hisat2-build 01-augustus/genome.fa index/chi_masked hisat2 -p 20 -x index/chi_masked -1 rna-seq/leaf_ox_r1_1.fastq.gz -2 rna-seq/leaf_ox_r1_2.fastq.gz | samtools sort -@ 10 > 02-barker/leaf_ox_r1.bam & isat2 -p 20 -x index/chi_masked -1 rna-seq/ox_flower9_rep1_1.fastq.gz -2 rna-seq/ox_flower9_rep1_2.fastq.gz | samtools sort -@ 10 > 02-barker/ox_flower9.bam & hisat2 -p 20 -x index/chi_masked -1 rna-seq/ox_flower16_rep1_1.fastq.gz -2 rna-seq/ox_flower16_rep1_2.fastq.gz | samtools sort -@ 10 > 02-barker/ox_flower16.bam & 然后,以未屏蔽重复序列的参考序列和BAM文件作为输入,让BRAKER2(安装会稍显麻烦,因为依赖许多软件)进行预测。 braker.pl --gff3 --cores 50 --species=carhr --genome=chi_unmasked.fa --bam=02-barker/leaf_ox_r1.bam,02-barker/ox_flower16.bam,02-barker/ox_flower9.bam # --gff3: 输出GFF3格式 # --genome: 基因组序列 # --bam: 比对后的BAM文件,允许多个 # --cores: 处理核心数 最后会得到如下输出文件 hintsfile.gff: 从RNA-seq比对结果的BAM文件中提取,其中内含子用于训练GeneMark-EX, 使用所有特征训练AUGUSTUS GeneMark-ET/genemark.gtf: GeneMark-EX根据RNA-seq数据训练后预测的基因 augustus.hints.gff: AUGUSTUS输出文件 将augustus.hints.gff3和文章的注释文件(carhr38.gtf)比较,见下图: comparation 其实不难发现,在不考虑UTR区域情况下,两者的差别其实更多表现是基因数目上,其实也就是利用转录组数据推测结构的问题所在,没有覆盖的区域到底是真的没有基因,还是有基因结构只不过所用组织没有表达,或者说那个区域其实是假基因?此外,如果基因间隔区域很短,有时候还会错误地把两个不同的基因预测为一个基因。因此,应该注重RNA-seq数据在剪切位点识别和外显子边界确定的优势。 03-同源预测基因结构 同源预测(homology prediction)利用近缘物种已知基因进行序列比对,找到同源序列。然后在同源序列的基础上,根据基因信号如剪切信号、基因起始和终止密码子对基因结构进行预测,如下示意图: 相对于从头预测的“大海捞针”,同源预测相当于先用一块磁铁在基因组大海中缩小了可能区域,然后从可能区域中鉴定基因结构。在10年之前,当时RNA-seq还没有普及, 只有少部分物种才有EST序列和cDNA序列的情况下,这的确是一个比较好的策略,那么问题来了,现在还需要进行这一步吗,如果需要是出于那种角度考虑呢? 在同源预测上,目前看到的大部分基因组文章都是基于TBLASTN + GeneWise,这可能是因为大部分基因组文章都是国内做的,这些注释自然而言用的就是公司的流程,然后目前国内的公司大多数又和某一家公司有一些关系。不过最近的3010水稻泛基因组用的是MAKER, 感谢部分提到这部分工作是由M. Roa(Philippine Genome Center Core Facilities for Bioinformatics, Department of Science)做的,算是一股清流吧。当然我在看Cardamine hirsuta基因组注释文章,发现它们同源注释部分用的是GenomeThreader, 该工具在本篇文章成文时的3月之前又更新了。 GeneWise的网站说它目前由Ewan Birney维护,只不过不继续开发了,因为Guy Slater开发Exonerate解决了GeneWise存在的很多问题,并且速度快了1000倍。考虑到目前只有GeneWise能利用HMM根据蛋白找DNA,而且ENSEMBL的注释流程也有一些核心模块用到了它,所以作者依旧在缓慢的开发这个工具(自2.4.1已经10多年没有更新了),当然这个工具也是非常的慢。尽管这一步不会用到GeneWise作为我们的同源注释选项,但是我们可以尝试用GeneWise手工注释一个基因,主要步骤如下 第一步: 使用BLASTX,根据dna序列搜索到蛋白序列,只需要第一个最佳比对结果 第二步: 选择最佳比对的氨基酸序列 第三步: 将dna序列前后延长2kb,与氨基酸序列一并传入给genewise进行同源预测 提取前5K序列,然后选择在TAIR上用BLASTX进行比对 seqkit faidx chi_unmasked.fa Chr1:1-5000 > chr1_5k.fa 让我们跳过这个尴尬的环节,毕竟很可能是我不太熟练使用工作所致。这里说点我的看法,除非你真的没有转录组数据,必须要用到同源物种的蛋白进行预测,或者你手动处理几个基因,否则不建议使用这个工具,因为你可能连安装都搞不定。 让我们用GenomeThreader基于上面的DNA序列和氨基酸序列进行同源基因结构预测吧 gth -genomic chr1_5k.fa -protein cer.fa -intermediate -gff3out # 其中cer.fa就是AT1G02205.2的氨基酸序列 结果一致,并且从RNA-seq的覆盖情况也符合预期 Chr1 gth exon 1027 1197 Parent=gene1 Chr1 MIPS_CARH_v3.8 exon 1027 1197 Chr1 gth exon 1275 1448 Parent=gene1 Chr1 MIPS_CARH_v3.8 exon 1275 1448 Chr1 gth exon 1541 1662 Parent=gene1 Chr1 MIPS_CARH_v3.8 exon 1555 1662 Chr1 gth exon 1807 2007 Parent=gene1 Chr1 MIPS_CARH_v3.8 exon 1807 2007 Chr1 gth exon 2085 2192 Parent=gene1 Chr1 MIPS_CARH_v3.8 exon 2085 2192 Chr1 gth exon 2294 2669 Parent=gene1 Chr1 MIPS_CARH_v3.8 exon 2294 2669 Chr1 gth exon 3636 3855 Parent=gene1 Chr1 MIPS_CARH_v3.8 exon 3636 3855 Chr1 gth exon 3971 4203 Parent=gene1 Chr1 MIPS_CARH_v3.8 exon 3971 4203 Chr1 gth exon 4325 4548 Parent=gene1 Chr1 MIPS_CARH_v3.8 exon 4325 4548 Chr1 gth exon 4676 4735 Parent=gene1 Chr1 MIPS_CARH_v3.8 exon 4676 4735 全基因组范围预测流程如下: 准备cDNA和或protein序列:在https://phytozome.jgi.doe.gov/p下载靠谱的物种的蛋白质序列,如 Arabidopsis thaliana, Oryza sativa, Brassica rapa, 查找文献寻找目前该物种的已有EST/cDNA序列,或者RNA-seq从头组装转录组。这里仅考虑用同源物种的蛋白序列进行比对分析,转录组从头组装数据用于PASA整体比对到参考基因组和更新已有的基因解雇。 分别测试下不同物种的同源注释结果 #run seperately gth -species arabidopsis -translationtable 1 -gff3 -intermediate -protein ~/db/protein_db/Athaliana_167_TAIR10.protein.fa.gz -genomic chi_unmasked.fa -o 03-genomethreader/Athaliana.gff3 & gth -species arabidopsis -translationtable 1 -gff3 -intermediate -protein ~/db/protein_db/BrapaFPsc_277_v1.3.protein.fa.gz -genomic chi_unmasked.fa -o 03-genomethreader/Brapa.gff3 & gth -species arabidopsis -translationtable 1 -gff3 -intermediate -protein ~/db/protein_db/Osativa_323_v7.0.protein.fa.gz -genomic chi_unmasked.fa -o 03-genomethreader/Osativa.gff3 & 在定性角度上来看,同源注释的结果和从头预测的没啥差别, 其中B. rapa和A. thaliana和C. hirsuta都属于十字花科,而O. sativa是禾本科, 所以前两者预测的效果好。 IGV展示 当然实际的同源注释流程中不能是单个物种分别预测,应该是将所有的蛋白序列进行合并,然后用BLASTX找到最优的联配,之后用GenomeThreader进行预测。PASA流程提到的UniRef90作为同源注释的搜索数据库可能是更好的选择,由于UniRef优先选择哪些人工审查、注释质量高、来源于模式动植物的蛋白,所以可靠性相对于直接使用同源物中可能更高。 BLASTX + GenomeThreader的代码探索中 04-RNA-seq的两种使用策略 对于RNA-seq数据,有两种使用策略,一种是使用HISAT2 + StringTie先比对再组装, 一种是从头组装,然后使用PASA将转录本比对到基因组上。 基于HISAT2 + StringTie 首先,使用HISAT2将RNA-seq数据比对到参考基因组, 这一步和之前相似,但是要增加一个参数--dta,使得StingTie能更好的利用双端信息 hisat2-build 01-augustus/genome.fa index/chi_masked hisat2 --dta -p 20 -x index/chi_masked -1 rna-seq/leaf_ox_r1_1.fastq.gz -2 rna-seq/leaf_ox_r1_2.fastq.gz | samtools sort -@ 10 > rna-seq/leaf_ox_r1.bam & hisat2 --dta -p 20 -x index/chi_masked -1 rna-seq/ox_flower9_rep1_1.fastq.gz -2 rna-seq/ox_flower9_rep1_2.fastq.gz | samtools sort -@ 10 > rna-seq/ox_flower9.bam & hisat2 --dta -p 20 -x index/chi_masked -1 rna-seq/ox_flower16_rep1_1.fastq.gz -2 rna-seq/ox_flower16_rep1_2.fastq.gz | samtools sort -@ 10 > rna-seq/ox_flower16.bam & samtools merge -@ 10 rna-seq/merged.bam rna-seq/leaf_ox_r1.bam rna-seq/ox_flower9.bam rna-seq/ox_flower16.bam 然后用StringTie进行转录本预测 stringtie -p 10 -o rna-seq/merged.gtf rna-seq/merged.bam 对于后续的EvidenceModeler而言,它不需要UTR信息,只需要编码区CDS,需要用TransDecoder进行编码区预测 util/cufflinks_gtf_genome_to_cdna_fasta.pl merged.gtf input/chi_masked.fa > transcripts.fasta util/cufflinks_gtf_to_alignment_gff3.pl merged.gtf > transcripts.gff3 TransDecoder.LongOrfs -t transcripts.fasta TransDecoder.Predict -t transcripts.fasta util/cdna_alignment_orf_to_genome_orf.pl \ transcripts.fasta.transdecoder.gff3 \ transcripts.gff3 \ transcripts.fasta > transcripts.fasta.transdecoder.genome.gff3 最后结果transcripts.fasta.transdecoder.gff3用于提供给EvidenceModeler 基于PASA 在多年以前,那个基因组组装还没有白菜价,只有几个模式物种基因组的时代,对于一个未测序的基因组,研究者如果要研究某一个基因的功能,大多会通过同源物种相似基因设计PCR引物,然后去扩增cDNA. 如果是一个已知基因组的物种,如果要大规模识别基因, 研究者通常会使用EST(expressed sequence tags)序列。 相对于基于算法的从头预测,cDNA和EST序列更能够真实的反应出一个基因的真实结构,如可变剪切、UTR和Poly-A位点。PASA(Progam to Assemble Spliced Alignments)流程最早用于拟南芥基因组注释,最初的设计是通过将全长(full-length)cDNA和EST比对到参考基因组上,去发现和更新基因组注释。其中FL-cDNA和EST序列对最后结果的权重不同。 PASA流程示意 这是以前的故事,现在的故事是二代转录组以及一些三代转录组数据,那么如何处理这些数据呢?我认为三代转录组相对于过去的FL-cDNA,而二代转录组数据经过拼接后可以看作是更长的EST序列。由于目前最普及的还是普通的mRNA-seq, 也就只介绍这部分流程。 考虑到我还没有研究过三代的全长转录组,分析过数据,这里的思考极有可能出错,后续可能会修改这一部分思考。 转录组组装使用Trinity(conda安装) cd rna-seq Trinity --seqType fq --CPU 50 --max_memory 64G --left leaf_ox_r1_1.fastq.gz,ox_flower16_rep1_1.fastq.gz,ox_flower9_rep1_1.fastq.gz --right leaf_ox_r1_2.fastq.gz,ox_flower16_rep1_2.fastq.gz,ox_flower9_rep1_2.fastq.gz & PASA是由30多个命令组成的流程,相关命令位于PASApipeline/scripts,为了适应不同的分析,有些参数需要通过修改配置文件更改, cp ~/opt/biosoft/PASApipeline/pasa_conf/pasa.alignAssembly.Template.txt alignAssembly.config # 修改如下内容 DATABASE=database.sqlite validate_alignments_in_db.dbi:--MIN_PERCENT_ALIGNED=80 validate_alignments_in_db.dbi:--MIN_AVG_PER_ID=80 上述几行配置文件表明SQLite3数据库的名字,设置了scripts/validate_alignments_in_db.dbi的几个参数, 表示联配程度和相似程度。后续以Trinity组装结果和参考基因组作为输入,运行程序: ~/opt/biosoft/PASApipeline/scripts/Launch_PASA_pipeline.pl -c alignAssembly.config -C -R -g ../chi_unmasked.fa -t ../rna-seq/trinity_out_dir/Trinity.fasta --ALIGNERS blat,gmap 最后结果如下: database.sqlite.pasa_assemblies_described.txt database.sqlite.pasa_assemblies.gff3 database.sqlite.pasa_assemblies.gtf database.sqlite.pasa_assemblies.bed 其中gff3格式用于后续的分析。 目前的一些想法, 将从头组装的转录本比对到参考基因组上很大依赖组装结果,所以和EST序列和cDNA相比,质量上还有一点差距。 05-整合预测结果 从头预测,同源注释和转录组整合都会得到一个预测结果,相当于收集了大量证据,下一步就是通过这些证据定义出更加可靠的基因结构,这一步可以通过人工排查,也可以使用EVidenceModeler(EVM). EVM只接受三类输入文件: gene_prediction.gff3: 标准的GFF3格式,必须要有gene, mRNA, exon, CDS这些特征,用EVidenceModeler-1.1.1/EvmUtils/gff3_gene_prediction_file_validator.pl验证 protein_alignments.gff3: 标准的GFF3格式,第9列要有ID信和和target信息, 标明是比对结果 transcript_alignments.gff3:标准的GFF3格式,第9列要有ID信和和target信息,标明是比对结果 EVM对gene_prediction.gff3有特殊的要求,就是GFF文件需要反映出一个基因的结构,gene->(mRNA -> (exon->cds(?))(+))(+), 表示一个基因可以有多个mRNA,即基因的可变剪接, 一个mRNA都可以由一个或者多个exon(外显子), 外显子可以是非翻译区(UTR),也可以是编码区(CDS). 而GlimmerHMM, SNAP等 这三类根据人为经验来确定其可信度,从直觉上就是用PASA根据mRNA得到的结果高于从头预测。 第一步:创建权重文件,第一列是来源类型(ABINITIO_PREDICTION, PROTEIN, TRANSCRIPT), 第二列对应着GFF3文件的第二列,第三列则是权重.我这里用了三个来源的数据。 mkdir 05-EVM && cd 05-EVM #vim weights.txt ABINITIO_PREDICTION augustus 4 TRANSCRIPT assembler-database.sqlite 7 OTHER_PREDICTION transdecoder 8 我觉得根据基因组引导组装的ORF的可信度高于组装后比对,所以得分和PASA差不多一样高。从头预测权重一般都是1,但是BRAKER可信度稍微高一点,可以在2~5之间。 第二步:分割原始数据, 用于后续并行. 为了降低内存消耗,--segmentsSize设置的大小需要少于1Mb(这里是100k), --overlapSize的不能太小,如果数学好,可用设置成基因平均长度加上2个标准差,数学不好,就设置成10K吧 cat transcripts.fasta.transdecoder.genome.gff3 ../braker/carhr/augustus.hints.gff3 > gene_predictions.gff3 ln ../04-align-transcript/database.sqlite.pasa_assemblies.gff3 transcript_alignments.gff3 ~/opt/biosoft/EVidenceModeler-1.1.1/EvmUtils/partition_EVM_inputs.pl --genome ../chi_unmasked.fa --gene_predictions gene_predictions.gff3 --transcript_alignments transcript_alignments.gff3 --segmentSize 100000 --overlapSize 10000 --partition_listing partitions_list.out 第三步:创建并行运算命令并且执行 ~/opt/biosoft/EVidenceModeler-1.1.1/EvmUtils/write_EVM_commands.pl --genome ../chi_unmasked.fa --weights `pwd`/weights.txt \ --gene_predictions gene_predictions.gff3 \ --transcript_alignments transcript_alignments.gff3 \ --output_file_name evm.out --partitions partitions_list.out > commands.list parallel --jobs 10 < commands.list 第四步:合并并行结果 ~/opt/biosoft/EVidenceModeler-1.1.1/EvmUtils/recombine_EVM_partial_outputs.pl --partitions partitions_list.out --output_file_name evm.out 第五步:结果转换成GFF3 ~/opt/biosoft/EVidenceModeler-1.1.1/EvmUtils/convert_EVM_outputs_to_GFF3.pl --partitions partitions_list.out --output evm.out --genome ../chi_unmasked.fa find . -regex ".*evm.out.gff3" -exec cat {} \; | bedtools sort -i - > EVM.all.gff 当前权重设置下,EVM的结果更加严格,需要按照实际情况调整,增加其他证据。 06-可选步骤 注释过滤:对于初步预测得到的基因,还可以稍微优化一下,例如剔除编码少于50个AA的预测结果,将转座子单独放到一个文件中(软件有TransposonPSI)。 这里基于gffread先根据注释信息提取所有的CDS序列,过滤出长度不足50AA的序列,基于这些序列过滤原来的的注释 gffread EVM.all.gff -g input/genome.fa -y tr_cds.fa bioawk -c fastx '$seq < 50 {print $comment}' tr_cds.fa | cut -d '=' -f 2 > short_aa_gene_list.txt grep -v -w -f short_aa_gene_list.txt EvM.all.gff > filter.gff 使用PASA更新EVM结果:EVM结果不包括UTR区域和可变剪切的注释信息,可以使用PASA进行更新。然而这部分已经无法逃避MySQL, 服务器上并没有MySQL的权限,我需要学习Perl脚本进行修改。因此基因结构注释到此先放一放。 07-基因编号 对每个基因实现编号,形如ABCD000010的效果,方便后续分析。如下代码是基于EVM.all.gff,使用方法为python gffrename.py EVM_output.gff prefix > renamed.gff. #!/usr/bin/env python3 import re import sys if len(sys.argv) < 3: sys.exit() gff = open(sys.argv[1]) prf = sys.argv[2] count = 0 mRNA = 0 cds = 0 exon = 0 print("##gff-version 3.2.1") for line in gff: if not line.startswith("\n"): records = line.split("\t") records[1] = "." if re.search(r"\tgene\t", line): count = count + 10 mRNA = 0 gene_id = prf + str(count).zfill(6) records[8] = "ID={}".format(gene_id) elif re.search(r"\tmRNA\t", line): cds = 0 exon = 0 mRNA = mRNA + 1 mRNA_id = gene_id + "." + str(mRNA) records[8] = "ID={};Parent={}".format(mRNA_id, gene_id) elif re.search(r"\texon\t", line): exon = exon + 1 exon_id = mRNA_id + "_exon_" + str(exon) records[8] = "ID={};Parent={}".format(exon_id, mRNA_id) elif re.search(r"\tCDS\t", line): cds = cds + 1 cds_id = mRNA_id + "_cds_" + str(cds) records[8] = "ID={};Parent={}".format(cds_id, mRNA_id) else: continue print("\t".join(records)) gff.close() 如果有转录组数据,没必须要使用太多的从头预测工具,braker2 加 GlimmerHMM可能就够用了, 更多是使用PASA和StringTie利用好转录组数据进行注释。 基因功能注释 基因功能的注释依赖于上一步的基因结构预测,根据预测结果从基因组上提取翻译后的 蛋白序列 和主流的数据库进行比对,完成功能注释。常用数据库一共有以几种: Nr:NCBI官方非冗余蛋白数据库,包括PDB, Swiss-Prot, PIR, PRF; 如果要用DNA序列,就是nt库 Pfam: 蛋白结构域注释的分类系统 Swiss-Prot: 高质量的蛋白数据库,蛋白序列得到实验的验证 KEGG: 代谢通路注释数据库. GO: 基因本体论注释数据库 除了以上几个比较通用的数据库外,其实还有很多小众数据库,应该根据课题研究和背景进行选择。注意,数据库本身并不能进行注释,你只是通过序列相似性进行搜索,而返回的结果你称之为注释。因此数据库和搜索工具要进行区分,所以你需要单独下载数据库和搜索工具,或者是同时下载包含数据库和搜索工具的安装包。 注意,后续分析中一定要保证你的蛋白序列中不能有代表氨基酸字符以外的字符,比如说有些软件会把最后一个终止密码子翻译成"."或者"*" BLASTP 这一部分用到的数据库都是用BLASTP进行检索,基本都是四步发:下载数据库,构建BLASTP索引,数据库检索,结果整理。其中结果整理需要根据BLASTP的输出格式调整。 Nr的NCBI收集的最全的蛋白序列数据库,但是无论是用NCBI的BLAST还是用速度比较快DIAMOND对nr进行搜索,其实都没有利用好物种本身的信息。因此在RefSeq上下载对应物种的蛋白序列, 用BLASTP进行注释即可。 # download wget -4 -nd -np -r 1 -A *.faa.gz ftp://ftp.ncbi.nlm.nih.gov/refseq/release/plant/ mkdir -p ~/db/RefSeq zcat *.gz > ~/db/RefSeq/plant.protein.faa # build index ~/opt/biosoft/ncbi-blast-2.7.1+/bin/makeblastdb -in plant.protein.faa -dbtype prot -parse_seqids -title RefSeq_plant -out plant # search ~/opt/biosoft/ncbi-blast-2.7.1+/bin/blastp -query protein.fa -out RefSeq_plant_blastp.xml -db ~/db/RefSeq/uniprot_sprot.fasta -evalue 1e-5 -outfmt 5 -num_threads 50 & Swiss-Prot里收集了目前可信度最高的蛋白序列,一共有55w条记录,数据量比较小, # download wget -4 -q ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz gzip -d uniprot_sprot.fasta.gz # builid index ~/opt/biosoft/ncbi-blast-2.7.1+/bin/makeblastdb -in uniprot_sprot.fasta -dbtype prot -title swiss_prot -parse_seqids # search ~/opt/biosoft/ncbi-blast-2.7.1+/bin/blastp -query protein.fa -out swiss_prot.xml -db ~/db/swiss_prot/uniprot_sprot.fasta -evalue 1e-5 -outfmt 5 -num_threads 50 & 关于结果整理,已经有很多人写了脚本,比如说我搜索BLAST XML CSV,就找到了https://github.com/Sunhh/NGS_data_processing/blob/master/annot_tools/blast_xml_parse.py, 所以就不过多介绍。 InterProScan 下面介绍的工具是InterProScan, 从它的9G的体量就可以感受它的强大之处,一次运行同时实现多个信息注释。 InterPro注释 Pfam数据库注释(可以通过hmmscan搜索pfam数据库完成) GO注释(可以基于NR和Pfam等数据库,然后BLAST2GO完成,) Reactome通路注释,不同于KEGG ./interproscan-5.29-68.0/interproscan.sh -appl Pfam -f TSV -i sample.fa -cpu 50 -b sample -goterms -iprlookup -pa -appl告诉软件要执行哪些数据分析,勾选的越多,分析速度越慢,Pfam就行。 KEGG数据库目前本地版收费,在线版收费,所以只能将蛋白序列在KEGG服务器上运行。因此你需要在http://www.genome.jp/tools/kaas/选择合适的工具进行后续的分析。我上传的50M大小蛋白序列,在KEGG服务器上只需要运行8个小时,也就是晚上提交任务,白天回来干活。 PASA:真核生物基因的转录本可变剪切自动化注释项目,需要提供物种的EST或RNA-seq数据 MAKER BRAKER1: 使用GeneMark-ET和AUGUSTUS基于RNA-Seq注释基因结构 EuGene JBrowse/GBrowse 参考文献和推荐阅读: NCBI真核生物基因组注释流程https://www.ncbi.nlm.nih.gov/genome/annotation_euk/process/ 真核基因组注释入门: "A beginner’s guide to eukaryotic genome annotation" 二代测序注释流程:Comparative Gene Finding: "Annotation Pipelines for Next-Generation Sequencing Projects" 基因组转录组注释策略: "Plant genome and transcriptome annotations: from misconceptions to simple solution" 重复序列综述: "Repetitive DNA and next-generation sequencing: computational challenges and solutions" MAKER2教程: http://weatherby.genetics.utah.edu/MAKER/wiki/index.php/MAKER_Tutorial_for_WGS_Assembly_and_Annotation_Winter_School_2018 《生物信息学》 樊龙江: 第1-5章: 基因预测与功能注释 《NGS生物信息分析》 陈连福: 真核生物基因组基因注释 JGS流程: https://genome.jgi.doe.gov/programs/fungi/FungalGenomeAnnotationSOP.pdf # Cardamine hirsutat基因组数据 mkdir chi_annotation && cd chi_annotation wget http://chi.mpipz.mpg.de/download/sequences/chi_v1.fa cat chi_v1.fa | tr 'atcg' 'ATCG' > chi_unmasked.fa # Cardamine hirsutat转录组数据 wget -4 -q -A '*.fastq.gz' -np -nd -r 2 http://chi.mpipz.mpg.de/download/fruit_rnaseq/cardamine_hirsuta/ & wget -4 -q -A '*.fastq.gz' -np -nd -r 2 http://chi.mpipz.mpg.de/download/leaf_rnaseq/cardamine_hirsuta/ & RepeatMasker: 用于注释基因组的重复区,需要安装RMBlast, TRF,以及在http://www.girinst.org注册以下载Repbase 安装RepeatMasker cd ~/src wget http://tandem.bu.edu/trf/downloadstrf409.linux64 mv trf409.linux64 ~/opt/bin/trf chmod a+x ~/opt/bin/trf # RMBlast cd ~/src wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.6.0/ncbi-blast-2.6.0+-src.tar.gz wget http://www.repeatmasker.org/isb-2.6.0+-changes-vers2.patch.gz tar xf ncbi-blast-2.6.0+-src gunzip isb-2.6.0+-changes-vers2.patch.gz cd ncbi-blast-2.6.0+-src patch -p1 < ../isb-2.6.0+-changes-vers2.patch cd c++ ./configure --with-mt --prefix=~/opt/biosoft/rmblast --without-debug && make && make install # RepeatMasker cd ~/src wget http://repeatmasker.org/RepeatMasker-open-4-0-7.tar.gz tar xf RepeatMasker-open-4-0-7.tar.gz mv RepeatMasker ~/opt/biosoft/ cd ~/opt/biosoft/RepeatMasker ## 解压repbase数据到Libraries下 ## 配置RepatMasker perl ./configure 在上面的基础上安装RepeatModel # RECON cd ~/src wget -4 http://repeatmasker.org/RepeatModeler/RECON-1.08.tar.gz tar xf RECON-1.08.tar.gz cd RECON-1.08/src make && make install cd ~/src mv RECON-1.08 ~/opt/biosoft # nesg cd ~/src mkdir nesg && cd nesg wget -4 ftp://ftp.ncbi.nih.gov/pub/seg/nseg/* mv nmerge nseg ~/opt/bin/ # RepeatScout http://www.repeatmasker.org/RepeatScout-1.0.5.tar.gz # RepeatModel wget -4 http://repeatmasker.org/RepeatModeler/RepeatModeler-open-1.0.11.tar.gz tar xf RepeatModeler-open-1.0.11.tar.gz mv RepeatModeler-open-1.0.11 ~/opt/biosoft/ cd ~/opt/biosoft/RepeatModeler-open-1.0.11 perl ./configure export PATH=~/opt/biosoft/maker:$PATH BLAST,BLAST有两个版本可供选择, WuBLAST或者NCBI-BLAST,我个人倾向于NCBI-BLAST,并且推荐使用编译后二进制版本,因为编译实在是太花时间了 cd ~/src wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.7.1+-x64-linux.tar.gz tar xf ncbi-blast-2.7.1+-x64-linux.tar.gz -C ~/opt/biosoft # 环境变量 export PATH=~/opt/biosoft/ncbi-blast-2.7.1+/bin:$PATH # 用于后续的BRAKER2 conda create -n annotation blast=2.2.31 AUGUSTUS: 可以说是最好的预测软件,使用conda安装 source activate annotation conda install augustus=3.3 GeneMark-ES/ET则是唯一一款支持无监督训练模型, 软件下载需要登记 cd ~/src wget http://topaz.gatech.edu/GeneMark/tmp/GMtool_Qg87n/gm_et_linux_64.tar.gz tar xf gm_et_linux_64.tar.gz mv gm_et_linux_64/gmes_petap/ ~/opt/biosoft wget http://topaz.gatech.edu/GeneMark/tmp/GMtool_Qg87n/gm_key_64.gz gzip -dc gm_key_64.gz > ~/.gm_key cpan YAML Hash::Merge Logger::Simple Parallel::ForkManager echo "export PATH=$PATH:~/opt/biosoft/gmes_petap/" >> ~/.bashrc GlimmerHMM: cd ~/src wget -4 ftp://ccb.jhu.edu/pub/software/glimmerhmm/GlimmerHMM-3.0.4.tar.gz tar xf GlimmerHMM-3.0.4.tar.gz -C ~/opt/biosoft SNAP: 基因从头预测工具,在处理含有长内含子上的基因组上表现欠佳 cd ~/src git clone https://github.com/KorfLab/SNAP.git cd SNP cd .. mv SNAP ~/opt/biosoft # 环境变量 export Zoe=~/opt/biosoft/SNAP/Zoe export PATH=~/opt/biosoft/SNAP:$PATH Exnerate 2.2: 配对序列比对工具,提供二进制版本, 功能类似于GeneWise,能够将cDNA或蛋白以gao align的方式和基因组序列联配。 cd ~/src wget http://ftp.ebi.ac.uk/pub/software/vertebrategenomics/exonerate/exonerate-2.2.0-x86_64.tar.gz tar xf exonerate-2.2.0-x86_64.tar.gz mv exonerate-2.2.0-x86_64 ~/opt/biosoft/exonerate-2.2.0 # .bashrc添加环境变量 export PATH=~/opt/biosoft/exonerate-2.2.0:$PATH conda install -c bioconda exonerate GenomeThreader 1.70: 同源预测软件,1.7.0版本更新于2018年2月 wget -4 http://genomethreader.org/distributions/gth-1.7.0-Linux_x86_64-64bit.tar.gz tar xf gth-1.7.0-Linux_x86_64-64bit.tar.gz -C ~/opt/biosoft # 修改.bashrc增加如下行 export PATH=$PATH:$HOME/opt/biosoft/gth-1.7.0-Linux_x86_64-64bit/bin export BSSMDIR="$HOME/opt/biosoft/gth-1.7.0-Linux_x86_64-64bit/bin/bssm" export GTHATADIR="$HOME/opt/biosoft/gth-1.7.0-Linux_x86_64-64bit/bin/gthdata" BRAKER2: 依赖AUGUSTUS 3.3, GeneMark-EX 4.33, BAMTOOLS 2.5.1, NCBI BLAST+ 2.2.31+(可选 SAMTOOLS 1.74+, GenomeThreader 1.70) cpan File::Spec::Functions Module::Load::Conditional POSIX Scalar::Util::Numeric YAML File::Which Logger::Simple Parallel::ForkManager cd ~/src wget -4 http://exon.biology.gatech.edu/GeneMark/Braker/BRAKER2.tar.gz tar xf BRAKER2.tar.gz -C ~/opt/biosoft echo "export PATH=$PATH:$HOME/opt/biosoft/BRAKER_v2.1.0/" >> ~/.bashrc # 在~/.bashrc设置如下软件所在环境变量 export AUGUSTUS_CONFIG_PATH=$HOME/miniconda3/envs/annotation/config/ export AUGUSTUS_SCRIPTS_PATH=$HOME/miniconda3/envs/annotation/bin/ export BAMTOOLS_PATH=$HOME/miniconda3/envs/annotation/bin/ export GENEMARK_PATH=$HOME/opt/biosoft/gmes_petap/ export SAMTOOLS_PATH=$HOME/miniconda3/envs/annotation/bin/ export ALIGNMENT_TOOL_PATH=$HOME/opt/biosoft/gth-1.7.0-Linux_x86_64-64bit/bin/ TransDecoder 编码区域预测工具,需要预先安装NCBI-BLAST cpan URI::Escape cd ~/src wget -4 https://github.com/TransDecoder/TransDecoder/archive/TransDecoder-v5.3.0.zip unzip TransDecoder-v5.3.0.zip cd TransDecoder-v5.3.0 make test MARKER: 使用conda安装会特别的方便,最好新建环境 conda create -n marker marker PASA: 依赖于一个数据库(MySQL或SQLite), Perl模块(DBD::mysql或DBD::SQLite), GMAP, BLAT, Fasta3。由于MySQL在HPC集群中的表现不如SQLite,以及安装MySQL还需要各种管理员权限,于是就有人进行了修改,增加了feature/sqlite分支, 见Add support for SQLite cpan DB_File URI::Escape DBI DBD::SQLite # GMAP wget http://research-pub.gene.com/gmap/src/gmap-gsnap-2017-11-15.tar.gz tar xf gmap-gsnap-2017-11-15.tar.gz cd gmap-2017-11-15 ./configure --prefix=$HOME/opt/gmap make && make install # BLAT cd ~/opt/bin wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/blat/blat && chmod 755 ./blat # Fasta3 wget -4 http://faculty.virginia.edu/wrpearson/fasta/fasta36/fasta-36.3.8g.tar.gz && \ tar zxvf fasta-36.3.8g.tar.gz && \ cd ./fasta-36.3.8g/src && \ make -f ../make/Makefile.linux_sse2 all && \ cp ../bin/fasta36 ~/opt/bin # 以上程序需添加到环境变量中 # PASApipeline cd ~/opt/biosoft git clone https://github.com/PASApipeline/PASApipeline.git cd PASApipeline && \ git checkout feature/sqlite && \ git submodule init && git submodule update && \ EVidenceModeler: 整合不同来源的注释结果,找到可靠的基因结构 cd ~/src wget -4 https://github.com/EVidenceModeler/EVidenceModeler/archive/v1.1.1.tar.gz tar xf v1.1.1.tar.gz mv EVidenceModeler-1.1.1 ~/opt/biosoft/

maker 在基因组注释上,MAKER算是一个很强大的分析流程。能够识别重复序列,将EST和蛋白序列比对到基因组,进行从头预测,并在最后整合这三个结果保证结果的可靠性。此外,MAKER还可以不断训练,最初的输出结果可以继续用作输入训练基因预测的算法,从而获取更高质量的基因模型。 Maker的使用比较简单,在软件安装成后,会有一个"data"文件夹存放测试数据 ls ~/opt/biosoft/maker/data dpp_contig.fasta dpp_est.fasta dpp_protein.fasta hsap_contig.fasta hsap_est.fasta hsap_protein.fasta te_proteins.fasta 以"dpp"开头的数据集为例,protein表示是同源物种的蛋白序列,est是表达序列标签,存放的是片段化的cDNA序列,而contig则是需要被预测的基因组序列。 让我们新建一个文件夹,并将这些测试数据拷贝过来。 mkdir test01 ; cd test01 cp ~/opt/biosoft/maker/data/dpp* . 由于基因组注释设计到多个程序,多个步骤,每个步骤可能都有很多参数需要调整,因此就需要建立专门的配置文件用来告诉maker应该如何控制流程的运行。 如下步骤创建三个以ctl结尾的配置文件 ~/opt/biosoft/maker/bin/maker -CTL ls *.ctl maker_bopts.ctl maker_exe.ctl maker_opts.ctl maker_exe.ctl: 执行程序的路径 maker_bopt.ctl: BLAST和Exonerat的过滤参数 maker_opt.ctl: 其他信息,例如输入基因组文件 maker_exe.ctl和maker_bopt.ctl可以简单用less查看,可不做修改,maker_opt.ctl是主要调整的对象。 使用vim maker_opt.ctl修改如下内容 genome=dpp_contig.fasta est=dpp_est.fasta protein=dpp_protein.fasta est2genome=1 修改完之后多花几分钟看看每个参数的设置,尽管很枯燥,但是考虑这个工具你可能会反复多次使用,所以这点时间是一定要花的。 随后就可以在当前路径运行程序 ~/opt/biosoft/maker/bin/maker &> maker.log & 输出结果见"dpp_contig.maker.output", 重点是"dpp_contig_master_datastore_index.log"文件,由于maker会拆分数据集并行计算,因此该文件记录总体的运行情况,需要关注其中是否有"FAILED","RETRY","SKIPPED_SAMLL","DIED_SIPPED_PERMANET",因为这意味着有些数据出于某些原因没有运算。 最后,我们需要将并行运算的结果进行整合,导出GFF文件, 转录本序列和蛋白序列 ~/opt/biosoft/maker/bin/fasta_merge -d dpp_contig_master_datastore_index.log ~/opt/biosoft/maker/bin/gff3_merge -d dpp_contig_master_datastore_index.log 在该目录下就会出现, "dpp_contig.all.gff", "dpp_contig.all.maker.proteins.fasta","dpp_contig.all.maker.transcripts.fasta" 其中GFF文件就需要用IGV,JBrowse, Apollo下展示来检查下注释是否正确。 软件安装:MAKER可以免费用于学术用途,但是未经许可不可商用。目前有两个版本2018年5月4日更新的2.31.10和测试版3.01.02.出于稳定性考虑,安装前者。后续假设已经在http://yandell.topaz.genetics.utah.edu/cgi-bin/maker_license.cgi进行登记,并且下载了压缩包"maker-2.31.10.tgz" 先检查下自己的系统情况,看需要补充哪些库 tar xf maker-2.31.10.tgz cd maker/src perl Build.PL 这一步之后会罗列出后续需要运行的命令来完成安装 ./Build installdeps ./Build installexes ./Build install ./Build status

如何解决matplotlib运行出现的Invalid DISPLAY variable

最近在服务器上运行matplotlib相关的脚本时遇到了"Invalid DISPLAY variable"报错,从报错中就可以知道这是因为没有显示设备导致的报错。 解决方案: 方案一: ~/.config/matplotlib/matplotlibr,在里面添加backend : Agg 这个方案不一定有用,如果失效考虑下面两种 方案二: 更换后端 可以先设置后端,然后导入pyplot import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt 或者先导入pyplot,然后切换后端 import matplotlib.pyplot as plt plt.switch_backend('Agg') 之后可以以Pdf形式或者其他格式保存到硬盘上。 from matplotlib.backends.backend_pdf import PdfPages import matplotlib.pyplot as plt plt.switch_backend('Agg') pdf = PdfPages('cut_figure.pdf') #先创建一个pdf文件 plt.figure() pdf.savefig() #将图片保存在pdf文件中 plt.close() pdf.close() #这句必须有,否则程序结束pdf文件无法打开

$PATH=~/.lib:$PATH; export PATH 当时我的回答是,"shell就是这样子规定的呀"。 回答的同时,也突然间发现有些自己感觉很熟悉的概念,原来自己并没有那么清楚,因此这一篇讲讲shell的命令行替换。先说结论 shell会在命令执行前对命令行进行一些替换 shell替换有如下几种: 之前使用命令 指定的文本 ~username 用户的主目录 $, ${...} shell和环境变量 $((..)) `...`, $(...) 运行在子shell里命令的输出 *,?,[...],[^...] 文件系统中匹配的文件名 历史替换是以!开头的替换方式,以下面历史记录为例 $ !-2 # 执行倒数第二个命令,即ls 大括号替换: 它会将{...}里的内容展开为多个单词,可以快速创建有一定规律的文件. 下面这个命令就把"chap0{1..3}"替换成了chap01, chap02, chap03, 以及每个都还有一个html和text对应。 $ mkdir -p chap0{1..3}/{html,text} $ tree chap0* chap01 ├── html └── text chap02 ├── html └── text chap03 ├── html └── text 代字号代替: 我们经常会看到别人文章会写用vim ~/.bashrc修改家目录下的配置文件,其中~默认就会替代成自己家目录路径,可以用echo ~确认。 那么问题来了,如何我想快速到别人的家目录下,应该怎么操作。只要在~加上别人的用户名就行了。比如说我/home 下还有一个用户叫做abc, 那么查看它家目录下的内容就是 ls ~abc 注: ~a可以用tab补全成~abc哦 变量替换: shell会把${变量名}或者$变量名替换成变量所指代的具体字符,比如说我将abc指代为ls,那么shell就会将$abc解释成ls,然后执行ls abc=ls # Desktop bin biosoft blastdb miniconda3 ncbi 也就是$PATH=~/.lib:$PATH; export PATH报错的原因是,shell在执行命令前会把$PATH成原来PATH里的字符串,显然无法达到修改PATH的目的 算术替换: shell命令行支持整数型的数学运算,下面的运算都是可以的,但是就别拿100/2.5这种浮点运算为难shell了。 echo $((1+2)) echo $((1-2)) echo $((100*101)) echo $((100/50)) 命令替换:这个替换非常的实用,可以将shell命令的输入结果立刻作为输入,而不是额外创建一个变量命。有一个应用场景就是在的分析报告里加上完成时间点 touch reports.$(date +%d%b%Y).log 路径名替换:路径替换的语法就4种,*表示0或更多的任意字符,?表示一个任意字符,[...]表示括号内的字符之一,[^...]不包括括号内的字符 以上就是shell命令行替换的几种形式。当然为了再一次强调"shell会在命令执行前对命令行进行一些替换",下面举一个反面例子来说明下。 Linux的/etc目录下有很多以conf结尾的配置文件,我们可以用find命令快速的定位到它们。 find /etc -name *.conf 上面的命令看起来没啥毛病,但是只要多做一件事情,就会有报错哦 touch a.conf b.conf find /etc -name *.conf # 如下是报错 find: paths must precede expression: b.conf Usage: find [-H] [-L] [-P] [-Olevel] [-D help|tree|search|stat|rates|opt|exec|time] [path...] [expression] 你会不会好奇,明明是相同的命令,却又不同的境遇呢?让我解释下,在刚开始的时候,文件下面没有"a.conf","b.conf",尽管shell看到"*"会有一种进行通配的冲动,但是很可惜没有对象让它统配。后来我们创建了这两个文件,给shell找到通配的机会,于是实际执行的命令就成了 "find /etc -name a.conf b.conf"。 由于后面这两个是文件路径,不符合find的命令要求,就导致了报错。 其实报错还好,有些时候没有报错,程序运行得到错误的结果反而更惨 如何避免这种错误呢?我们就需要用到"避免\*这个元字符被shell解释。 除了双引号,避免shell进行替换的符号还有 反斜杠\ , 和单引号 '. 单引号和双引号的区别在于,单引号内部所有字符都是普通字符而已,而双引号里的美元符号$, 感叹号! 和反引号 ` 还能被shell解释

这下简书上的markdown完整了

最近用简书的markdown写作时,偶尔发现简书又支持几个markdown扩展语法。 第一个是上标, 实现方法是用两个^包围需要上标的字符,如hoptop^TM^, 效果为hoptopTM 第二个是下标, 实现方法是用两个~包围需要下标的字符,如hoptop~TM~, 效果就是hoptopTM 这两个语法在简书官方markdown推文-献给写作者的 Markdown 新手指南-中并没有提到, 算时最近写作的一些小惊喜吧。 更新于2018-9-3 简书终于支持数学符号了,激动呀 然而一些诡异的字符还需要慢慢加油哦, 比如说下面这几个 这个操作有一个缺点就是,如果我想从管道里面传入数据给Python的话,就会报错,因为原代码要求文件而不是标准输入。 这个问题可以通过Python的一个标准库: fileinput进行解决。 import fileinput for line in fileinput.input(): process(line) fileinput.intpu()会帮我们自动处理输入。如果sys.argv[1:]里有输入文件,它就会对里面所有的文件进行遍历,如果sys.argv为空,那么它就会从标准输入sys.stdin里读取输入,如果输入文件的文件名是"-", 同样地会从标准输入里读取输入。这样子就省去了我们自己写条件语句进行判断输入类型。

这是R数据科学的读书笔记之一,《R数据科学》是一本教你如何用R语言进行数据分析的书。即便我使用R语言快2年多了,但是读这本书还是受益颇多。 最早接触R语言的时候看的是《R语言实战》, 在第二章里,该书将R语言的数据结构分为6种,向量、矩阵、数组、数据框、因子和列表。当时的理解是,矩阵是二维的向量,数组是二维以上的向量,数据框是特殊性质的列表。 但是读完《R数据科学》的第15章:向量后,我发现原来R语言的数据结构原来可以只分为两类 原子向量: 包含6种类型,逻辑性、整型、双精度型、字符型、复数型和原始型 递归向量: 更常见的名字叫做列表 原子向量和递归向量的 唯一区别 就在于其中存放的值是否都是同种类型。 向量(vector), 矩阵(matrix)和数组(array)以及因子(factor)都只能存放一种数据类型,因此is.atomic的判断结果都是TRUE, 所以都是原子向量 数据库和列表可以包含不同类型的数据,所以用is.recursive的判断结果是TRUE,所以都是递归向量 此外,每个向量都有两个关键属性(properties),类型和长度, 分别用typeof()和length()进行查看。分别去用typeof()查看向量、矩阵、数组、因子、数据框和列表时,你会发现前面4个返回都是6种基本数据类型,而数据框和列表返回的都是"list". 我们还可以在向量上附加任意多的元数据(metadata),这些元数据称之为特征(attributes)。 附加不同的特性后就得到了扩展向量(augmented vectors), 其中名称、维度和类是三种特别重要的属性。 如果你去查看attribute和property的中文翻译时,你会发现两者都有一个释义叫做属性 从扩展向量的角度上看数据类型时,可以得到如下洞见 第一: 矩阵和数组相对于普通向量主要就多了一个dim属性,所以我们可以通过如下的操作来创建矩阵和数组 is.v.m.a <- function(x) {c(is.vector(x), is.matrix(x), is.array(x))} v <- c(1,2,3,4) is.v.m.a(v) # TRUE FALSE FALSE attr(v,'dim') <- c(2,2) is.v.m.a(v) # FALSE TRUE TRUE attr(v,'dim') <- c(1,2,2) is.v.m.a(v) # FALSE FALSE TRUE 注: 矩阵是特殊的数组。 第二:名称是一种额外属性, 对于向量是"names", 对于数组则是"dimnames[[x]]", x表示不同维度, 对于列表而言则是"names",对于数据框是"names"对于列名和"row.names"对于行名 v <- c(1,2,3,4) attr(v,'names') <- c('a','b','c','d') 第三:类(class)也是一种属性,类是面向对象编程的一个概念。在R语言中,我们会发现同一个函数居然可以用在不同的数据集,比如说print用在ggplot2的对象中,结果是输出图片,这种函数就称之为泛型函数。 methods(print)# 内容过多,不在这里展示 # 我们可以具体某个函数的代码 getS3method("print","data.frame") 关于泛型函数的更多知识会在后续的面向对象编程里介绍。 其他知识点 R语言的缺失值一般都标记为"NA", 因此在读取数据的时候默认也将文件中的"NA"当作缺失值,但是很有可能其他人会用"null"作为缺失值的标记,所以结果就会导致这一列全部被当做是字符串,影响后续的分析。 在向量取子集时,熟悉Python的人需要注意一点,Python中x=[1,2,3,4]; x[-1]表示选择最后一个元素,而在R语言里x= c(1,2,3,4); x[-1]表示删除第一个元素,即R用负整数取子集时会丢弃对应位置的元素。 [和[[在提取列表时,一定要注意,[[会使列表降低一个层次,而[会返回一个新的、更小的列表,也就是 l <- list(c(1,2,3)) l[1] # 返回列表 l[[1]] # 返回向量 为了更好理解这两者在列表中的差异,作者还提供了一个非常形象的例子,我用另一个例子来说明下: 我所就读的初中每个年级段大概有10个班级,每个班级的人数都不太一样。那么这里的一个年级段就是一个列表x,每个班级都是列表里元素。那么x[1]表示的是解散其他所有班级,只留下第一个班级组成年级段。而x[[1]]表示是第一个班级。x[[1]][1]表示的可能是第一个班级里的第一个学生。

使用Mikado进行基因结构注释

Mikado是基于Python3写的基因组结构注释工具,它主要做的是从多个转录组组装工具得到的转录本里挑选出最好的结果作为基因组的结构注释。此外,它还会基于同源蛋白比对结果对转录本打分。换句话说这个软件主要是根据转录组数据进行注释,没有 ab inito 预测。 软件安装比较方法,我们可以使用bioconda进行安装: conda create -n mikado mikado # 打开Python进行测试, 注意大小写 # import Mikado # Mikado.test() 使用Daijin准备Mikado所需文件 第一步: 准备输入 如下是下载参考序列和对应的GTF注释文件。 mkdir -p Reference cd Reference wget ftp://ftp.ensembl.org/pub/release-89/gtf/drosophila_melanogaster/Drosophila_melanogaster.BDGP6.89.gtf.gz wget ftp://ftp.ensembl.org/pub/release-89/fasta/drosophila_melanogaster/dna/Drosophila_melanogaster.BDGP6.dna.toplevel.fa.gz wget "http://www.uniprot.org/uniprot/?sort=score&desc=&compress=yes&query=taxonomy:diptera%20NOT%20taxonomy:%22Drosophila%20(fruit%20flies)%20[7215]%22%20AND%20taxonomy:%22Aedes%20aegypti%22&fil=&format=fasta&force=yes" -O Aedes_aegypti.fasta.gz gunzip *gz cd ../ 如下代码下载转录组数据 mkdir -p Reads cd Reads wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR166/003/ERR1662533/ERR1662533_1.fastq.gz wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR166/003/ERR1662533/ERR1662533_2.fastq.gz wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR166/004/ERR1662534/ERR1662534_1.fastq.gz wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR166/004/ERR1662534/ERR1662534_2.fastq.gz cd ../ 第二步:创建配置文件 使用daijin comfigure创建配置文件, 包括如下内容 配置文件名: -o OUT 每个任务的线程数: --threads N 物种名和参考序列:--name, --genome, --transcriptome 二代测序数据: --sample-sheet 比对软件: -al [{gsnap,star,hisat,tophat} [{gsnap,star,hisat,tophat} ...]] 组装工具: -as [{class,cufflinks,stringtie,trinity,scallop} [{class,cufflinks,stringtie,trinity,scallop} ...]] 输出文件夹: -od OUT_DIR 打分文件用于Mikado: --scoring 蛋白数据库: --prot-db 转录本的间距: --flank 集群任务投递工具: --scheduler 集群任务投递配置文件:-c CLUSTER_CONFIG daijin configure --scheduler "" \ --scoring dmelanogaster_scoring.yaml \ --copy-scoring dmelanogaster_scoring.yaml \ -m permissive --sample-sheet sample_sheet.tsv \ --flank 500 -i 50 26000 --threads 2 \ --genome Reference/Drosophila_melanogaster.BDGP6.dna.toplevel.fa \ -al hisat -as class stringtie -od Dmelanogaster -name Dmelanogaster \ -o daijin.yaml --prot-db Reference/Aedes_aegypti.fasta; 这里面的samples_sheet.tsv内容如下. 第一列和第二列是双端测序的read, 第三列是样本名, 第四列表示是否为链特异性建库, 包括非链特异性(fr-unstranded), 链特异性数据且第一个reads是正向链第二个reads是反向链(fr-firststrand ), 链特异性数据且第二个reads是正向链第一个reads是反向链(fr-secondstrand), 仅正向链(f)和仅反向链(r), 最后一列表示是否非为二代测序结果(False表示为二代测序) Reads/ERR1662533_1.fastq.gz Reads/ERR1662533_2.fastq.gz ERR1662533 fr-unstranded False Reads/ERR1662534_1.fastq.gz Reads/ERR1662534_2.fastq.gz ERR1662534 fr-unstranded False 第三步:运行 执行组装步骤。 daijin assemble --cores 20 -nd 运行时出现的问题和解决方案: 对于某些服务器而言,即便在参数将任务投递系统设置为空,程序依旧会投递,解决方案就是加上-nd 运行过程中会用到gnuplot进行绘图,如果报错,就找到对应行将其注释, 即下面的plot-bamstats部分。 rule bam_stats: input: bam=rules.bam_sort.output, idx=rules.bam_index.output output: ALIGN_DIR+"/output/{align_run}.sorted.bam.stats" params: load=loadPre(config, "samtools"), #plot_out=ALIGN_DIR+"/output/plots/{align_run}/{align_run}" threads: 1 message: "Using samtools to collected stats for: {input}" shell: "{params.load} samtools stats {input.bam} > {output}" #" && plot-bamstats -p {params.plot_out} {output}" 运行结束之后得到如下文件 Dmelanogaster/3-assemblies/output/class-0-hisat-ERR1662533-0.gtf Dmelanogaster/3-assemblies/output/class-0-hisat-ERR1662534-0.gtf Dmelanogaster/3-assemblies/output/stringtie-0-hisat-ERR1662533-0.gtf Dmelanogaster/3-assemblies/output/stringtie-0-hisat-ERR1662534-0.gtf 同时将总的统计结果存放在了"Dmelanogaster/3-assemblies/assembly.stats"下 第四步:运行Mikado 上一步提供了组装得到的GTF文件就可以作为Mikado的输入进行结构注释, 其中mikado要求的输入文件在dmelanogaster_scoring.yaml, 里面的内容 daijin mikado -nd Dmelanogaster/mikado.yaml 最后的结果在Dmelanogaster/5-mikado/pick/permissive/mikado-permissive.loci.gff3中 如果组装结果的GTF文件有一个为空那么就会报错,把这个组装软件在参数中删掉

不要轻易相信AnnotationHub的物种注释包

Bioconductor开发的物种注释包系列集合了一个物种不同来源的注释信息,能够根据基因ID对其进行多种来源的注释,比如说基因的别名,基因的功能等。 我之前也写过一篇文章用Bioconductor对基因组注释介绍如何使用AnnotationHub下载注释数据库, 使用select(), mapIds等函数进行注释操作。我自己写一个流程也用到了它给基因ID, 如AT1G14185, 注释别名和功能描述。 注释结果中会出现一些基因无法被注释, 比如说下面这些情况, 我一直认为只是这些基因没有得到比较好的研究, 即便这些基因能够在TAIR搜到Araport的注释, 我也认为那些注释只是同源注释没有多大意义。 AT1G13970 NA NA AT1G14120 NA NA AT1G14240 NA NA AT1G14600 NA NA 一开始得到的结果里没有多少个基因,所以缺少几个注释,通过手工去查找也行,但是目前差异表达分析动不动就给别人500多个基因,于是就有几十个甚至上百个未注释的基因,所以我想着要不自己更新拟南芥的物种包。 library(AnnotationHub) ah <- AnnotationHub() org <- ah[["AH57965"]] org#...# TAIRGENEURL: ftp://ftp.arabidopsis.org/home/tair/Genes/TAIR10_genome_release/TAIR10_functional_descriptions#... 通过上面的代码,我找到了基因功能描述的数据库来源文件,我下载了这个文件,并且拿用AnnotationHub注释不到的功能的一个基因,”AT1G14185”,进行测试 mapIds(org, "AT1G14185", "GENENAME", "TAIR") 这下就非常有趣了,在原始文件中能搜索到的基因用Bioconductor的物种注释包时却没有注释信息!为了搞清楚这个原因,我花了快半个下午的时间去折腾,终于被我找到了原因。 我分别用一个能被org.At.tair.db注释和一个不能被org.At.tair.db的注释去搜索原始文本. # 无注释 1 AT1G14185.1 2 protein_coding 3 Glucose-methanol-choline (GMC) oxidoreductase family protein 5 Glucose-methanol-choline (GMC) oxidoreductase family protein; FUNCTIONS IN: ... # 有注释 1 AT1G19610.1 2 protein_coding 3 Arabidopsis defensin-like protein 4 Predicted to encode a PR (pathogenesis-related) protein. ... 5 PDF1.4; FUNCTIONS IN: molecular_function unknown; 简单的比较之后,你差不多就知道了org.At.tair.db的在功能描述这一部分其实只用第一列和第四列(为了方便展示我转置了原始数据)。这就是非常让人意外了,为啥它不用第一列和第三列呢? 我于是又去看了其他几个基因,就差不多明白了,原始的文本特别的混乱,你除了能保证第一列和第二列有信息外,其他列你根本无法保证,因此最好的策略以第一列作为检索的关键字,其他列合并成一列才行,然而作者没有那么细致。 于是我就放弃了用org.At.tair.db注释基因功能描述和基因别名了,还是自己写一个Python脚本进行注释吧。 下面这个脚本只适用于bed格式的输入,且第四列为转录本ID,另外两个输入文件分别为”gene_aliases_20140331.txt”和”TAIR10_functional_descriptions_20140331.txt”, 用法为 python bed_anno.py to_anno.bed gene_aliases_20140331.txt TAIR10_functional_descriptions_20140331.txt > anno.xls import sysfrom collections import defaultdict bed_file = sys.argv[1] alias_file = sys.argv[2] func_file = sys.argv[3] alias_dict = defaultdict(list) func_dict = defaultdict(list) # read alias file for line in open(alias_file, 'r'): items = line.strip().split('\t') alias_dict[items[0]] = items[1:] # read function description file for line in open(func_file, 'r'): items = line.strip().split('\t') func_dict[items[0]] = items[1:] # annotation and output for line in open(bed_file, 'r'): transcript_id = line.strip().split("\t")[3] gene_id = transcript_id.split(".")[0] gene_alias = alias_dict[gene_id] if len(alias_dict[gene_id]) > 0 else [''] gene_func = func_dict[transcript_id] if len(func_dict[transcript_id]) > 0 else [''] gene_anno = '{}\t{}\t{}'.format(line.strip(), gene_alias[0], '\t'.join(gene_func)) print(gene_anno)

这是R数据科学的读书笔记之一,《R数据科学》是一本教你如何用R语言进行数据分析的书。即便我使用R语言快2年多了,但是读这本书还是受益颇多。 这一篇学习笔记对应第13章:使用magrittr进行管道操作。关于管道这个概念,我最早在Linux系统中接触,它是Unix系统设计哲学的体现,“组合小功能完成大任务”,比如说BWA比对后排序用管道的写法就是 bwa mem ref 1.fq 2.fq | samtools sort > align.bam 在R语言接触管道符号"%>%"是在学习dplyr包时候,那个时候我以为这个符号是 Hadley Wickham 创造出来的,其实是来源于Stefan Milton Bache开发的magrittr中。 在没有管道符号之前,如果我要对一个变量做一系列的分析的话,那么写法是下面这个样子 # 先创建100个随机数 nums <- rnorm(100) # 分成两列 nums_matrix <- matrix(nums, ncol = 2) # 分别求两列的均值 nums_mean <- Matrix::colMeans(nums_matrix) 这里面我写了很多中间变量,要多敲很多字,而且如果我要修改输入的话的100个随机数的话,我需要修改两处。当然可以进行函数嵌套. Matrix::colMeans(matrix(rnorm(100), ncol=2)) 但是这种写法不利于人的阅读,当我读到这个函数的时候,我需要先连续往大脑里塞进去两个函数后,才能抵达核心,然后再从里往外解析。 但是有了管道符号之后一切就不一样了,写法就是 rnorm(100) %>% matrix(ncol=2) %>% Matrix::colMeans() 你会发现从左往右阅读,代码读起来非常的流畅。 虽然管道看起来很美好,但是在如下的场景中就不太适合了, 操作步骤特别的多,比如说10个,那么你就需要用一些有意义的中间变量来存放中间结果,方便调试 多输入多输出。比如说A和B输入,输出C和D 操作步骤构成了一张复杂关系的有向图,比如说D结果依赖于B和C,而B和C依赖于A。 简单点说,就是类似于A > B > C > D 这种场景用管道比较好。 除了%>%这个好用的符号外,magrittr还提供了其他三个比较好用的符号,%$%,%<>%和%T>%。 上面都是常规操作,作为有一定基础的R语言使用者,更希望探索点这个符号的本质。 首先明确一点,在R语言中一切符号本质上都是函数,比如说"+"也是一个函数,常规用法都是1 + 2, 但是我们可以用函数的方式来写哦 `+`(4,5) 因此rnorm(100) %>% matrix(ncol=2)其实应该理解成 `%>%`(rnorm(100), matrix(ncol=2)) 那么我们就可以看看管道符号的源代码了 ?magrittr::`%>%` function (lhs, rhs) parent <- parent.frame() env <- new.env(parent = parent) chain_parts <- split_chain(match.call(), env = env) pipes <- chain_parts[["pipes"]] rhss <- chain_parts[["rhss"]] lhs <- chain_parts[["lhs"]] env[["_function_list"]] <- lapply(1:length(rhss), function(i) wrap_function(rhss[[i]], pipes[[i]], parent)) env[["_fseq"]] <- `class<-`(eval(quote(function(value) freduce(value, `_function_list`)), env, env), c("fseq", "function")) env[["freduce"]] <- freduce if (is_placeholder(lhs)) { env[["_fseq"]] else { env[["_lhs"]] <- eval(lhs, parent, parent) result <- withVisible(eval(quote(`_fseq`(`_lhs`)), env, env)) if (is_compound_pipe(pipes[[1L]])) { eval(call("<-", lhs, result[["value"]]), parent, parent) else { if (result[["visible"]]) result[["value"]] else invisible(result[["value"]]) 这个代码的核心在于如下两行 env[["_function_list"]] <- lapply(1:length(rhss), function(i) wrap_function(rhss[[i]], pipes[[i]], parent)) env[["_fseq"]] <- `class<-`(eval(quote(function(value) freduce(value, `_function_list`)), env, env), c("fseq", "function")) 这两行干的活其实是进行词法转换,也就是把我们之前的管道串联起来的部分转换成 my_pipe <- function(.){ . <- rnorm(.) . <- matrix(., ncol = 2) . <- Matrix::colMeans(.) my_pipe(100) 考虑到在管道里面用"+"."-"这些函数时用到`+`或许会有点诡异,于是magrittr给这些符号命名了对应的别名,如下 extract `[` extract2 `[[` inset `[<-` inset2 `[[<-` use_series `$` add `+` subtract `-` multiply_by `*` raise_to_power `^` multiply_by_matrix `%*%` divide_by `/` divide_by_int `%/%` mod `%%` is_in `%in%` and `&` or `|` equals `==` is_greater_than `>` is_weakly_greater_than `>=` is_less_than `<` is_weakly_less_than `<=` not (`n'est pas`) `!` set_colnames `colnames<-` set_rownames `rownames<-` set_names `names<-`

初步组装的杂合基因组如何去冗余

redundans的目标是辅助杂合基因组的组装,输入文件可以是组装的contig,测序文库以及额外的参考基因组,最后用于搭建出scaffold级别的纯合基因组组装结果。包括如下几个步骤: 从头组装: 它会调用Platanus、SSPACE3进行组装 去冗余: 从最初组装中去除冗余的序列 scaffolding: 利用双端测序将contig进行搭接 gap closing: 即填补scaffold中的N序列 对于我们三代组装的结果而言,我们只需要去冗余这一步即可。 这一步一定要保证你的电脑上装了ZLIB库,不然就需要去修改BWA和LAST的Makefile, 手动添加"CFLAGS"和"LDFLAGS", 你或许不行。 git clone --recursive https://github.com/lpryszcz/redundans.git cd redundans && bin/.compile.sh 结果输出"done"才算是成功.如果还需要作图,则需要安装 matplotlib numpy pip install matplotlib numpy 最好用下面这行命令测试下。 ./redundans.py -v -i test/*_?.fq.gz -f test/contigs.fa -o test/run1 软件的使用 这个软件就是在安装的时候让我折腾了下,使用倒是非常的方便,去冗余主要调整的参数就是相似度和重叠(overlap)度 默认相似度参数--identity 0.51,重叠比例是--overlap 0.80 越大越严格。 此外,如果你用-i参数提供了二代测序数据,redundans还会默认搭scaffold和补洞,但我只需要用到它的去冗余步骤, 另外的两步我不要,所以还要添加--noscaffolding和--nogapclosing跳过这两步。 ident=0.55 ovl=0.80 contig=contig.fa threads=10 redundans.py -v -f ${contig} -o ident_${ident}_ovl_${ovl} -t ${threads} \ --log ident_${ident}_ovl_${ovl}.log \ --identity ${ident} --overlap ${ovl} \ --noscaffolding --nogapclosing 上面代码运行时如果不小心中断了,加上--resume就能断点重跑了。 https://github.com/lpryszcz/redundans Redundans: an assembly pipeline for highly heterozygous genomes

如何使用GMAP/GSNAP进行转录组序列比对

GMAP最早用于讲EST/cDNA序列比对到参考基因组上,可以用于基因组结构注释。后来高通量测序时代,又开发了GSNAP支持高通量数据比对,这篇文章主要介绍GMAP,毕竟高通量转录组数据比对大家更喜欢用STAR, HISTA2等软件。 下面是我源码安装的代码 wget http://research-pub.gene.com/gmap/src/gmap-gsnap-2018-07-04.tar.gz tar xf gmap-gsnap-2018-07-04.tar.gz cd gmap-2018-07-04/ ./configure --prefix=$HOME/opt/biosoft/gmap make -j 20 如下步骤假设你有一个物种的基因组序列和对应的CDS序列,分别命名为"reference.fa"和"cds.fa" 第一步:构建GMAP/GSNAP索引数据库 GMAP/GSNAP对FASTA文件中每个记录下的序列的长度有一定限制, 每一条不能超过4G, 能应付的了大部分物种了。 构建索引分为两种情况考虑,第一种是一个fasta文件包含所有的序列 ~/opt/biosoft/gmap/bin/gmap_build -d reference reference.fa 第二种则是每个染色体的序列都单独存放在一个文件夹里,比如说你下载人类参考基因组序列解压后发现有N多个fasta文件, 然后你就想用其中几条染色体构建索引 ~/opt/biosoft/gmap/bin/gmap_build -d reference Chr1.fa Chr2.fa Chr3.fa ... 注: 这里的-d表示数据K库的名字,默认把索引存放在gmap安装路径下的share里,可以用-D更改.此外还有一个参数-k用于设置K-mer的长度, 默认是15, 理论上只有大于4GB基因组才会有两条一摸一样的15bp序列(当然是完全随机情况下)。 第二步:正式比对 建立完索引之后就可以将已有的CDS或者EST序列和参考基因组序列进行比较。 ~/opt/biosoft/gmap/bin/gmap -t 10 -d reference -f gff3_gene cds.fa > cds_gene.gff3 其中-t设置线程数, -d表示参考基因组数据库的名字, 都是常规参数。我比较感兴趣的参数是如何将序列输出成GFF格式. GMAP允许多种格式的输出,比如说-S只看联配的总体情况,而-A会显示每个比对上序列的联配情况, 还可以输出蛋白序列(-P)或者是genomic序列(-E). 但是做结构注释要的gff文件,参数就是-f gff3_gene, -f gff3_match_cdna, -f gff3_match_est。 要想对一个软件有更好的认识,最好还是看看他们文章是怎么说的。 GMAP: a genomic mapping and alignment program for mRNA and EST sequences Bioinformatics 2005 21:1859-1875 Abstract Full Text, Thomas D. Wu and Colin K. Watanabe Fast and SNP-tolerant detection of complex variants and splicing in short reads Bioinformatics 2010 26:873-881 AbstractFull Text, Thomas D. Wu and Serban Nacu

出现的问题 这个问题非常的有趣,因为我在两台服务器分别安装同一个包,其中一台没有任何问题,而另一台则失败,尽管操作系统都是“CentOS Linux release 7.4.1708 (Core)”。不过我知道,一定还有一些差异我没有发现,有可能是R的版本不同,有可能是安装R所用的GCC版本不同,但是这种差异就算知道了也不像去解决,我可不想为了一个R包重装系统。 那么如何解决这个问题呢?策略就是手动安装了。 让我们先下载这个R包并解压 wget https://mirrors.tuna.tsinghua.edu.cn/CRAN/src/contrib/expm_0.999-2.tar.gz tar xf expm_0.999-2.tar.gz 然后找到locale.h中报错行, 也就是从libintl.h中调用的dgettext报错了 #include <R.h> #ifdef ENABLE_NLS #include <libintl.h> #define _(String) dgettext ("expm", String) #else #define _(String) (String) #endif 这个libintl.h凭我本能的直觉我认为是应该是在/usr/include/下, 使用VIM打开并查找"LC_MESSAGES" /* We need LC_MESSAGES for `dgettext'. */ # include <locale.h> 这行代码告诉我们这个"LC_MESSAGES"在locale.h下, 但是我发现在locale.h里面是有定义的 /* These are the possibilities for the first argument to setlocale. The code assumes that the lowest LC_* symbol has the value zero. */ #define LC_CTYPE __LC_CTYPE #define LC_NUMERIC __LC_NUMERIC #define LC_TIME __LC_TIME #define LC_COLLATE __LC_COLLATE #define LC_MONETARY __LC_MONETARY #define LC_MESSAGES __LC_MESSAGES #define LC_ALL __LC_ALL #define LC_PAPER __LC_PAPER #define LC_NAME __LC_NAME #define LC_ADDRESS __LC_ADDRESS #define LC_TELEPHONE __LC_TELEPHONE #define LC_MEASUREMENT __LC_MEASUREMENT #define LC_IDENTIFICATION __LC_IDENTIFICATIO 于是我崩溃了。 突然间我灵机一闪,为啥不直接把这个定义添加到原来的locale.h里呢? 于是我在locale.h里增加了一行,改为 /* Localization */ #define LC_MESSAGES "en_US.UTF-8" #include <R.h> #ifdef ENABLE_NLS #include <libintl.h> #define _(String) dgettext ("expm", String) #else #define _(String) (String) #endif 最后用R CMD INSTALL安装 R CMD INSTALL expm/ 顺利安装!

这是一篇工具介绍贴,考虑这个工具是要钱的,那些动不动就说别人忘了初心的用户肯定认为我写的是软文,所以这些人就不要继续往下看了。 变异检测的软件目前虽然有很多,SAMtools/BCFtools, GATK, FreeBayes等,但是我看到的大部分文章都是用GATK UG/HC。GATK的速度是有目共睹的慢,不过我平时就分析几个重测序样品,基本上过个两天就能出结果,所以速度不是我的刚需。 如果想要追求速度的话,一种思路是可以将参考基因组进行分割,然后分别并行运算加速,或者搭建Spark环境,用GATK4的Spark模式。还有一种就是根据GATK的算法思想,用C/C++重新写软件。去年的时候我看到了一个软件叫做sentieon,用C/C++实现了GATK的算法,瞬间速度就上来了,这是一个商业公司的收费软件,目前国内用的比较少。 我原本以为这个速度已经够快的,直到我最近去demo了另一个软件,edico公司开发的DRAGEN,这个效率简直是丧心病狂。它从硬件和软件上同时进行加速 需要购买他们公司的硬件,128G内存,56线程,2T固态硬盘,以及一个FPGA芯片 比对这一步的算法基于而不是BWT转换(这就是为什么要128G内存) 比对之前要将索引加载到芯片中,所以每次只能比对一个任务 由于IO读写是非常大瓶颈,所以采用了固态硬盘 程序由C/C++开发,所以效率极高 为啥我要去demo这个工具呢,主要因为最近服务器资源紧缺(因为之前用的服务器要么是合租的,要么是蹭别人的),而老板又在催进度,而要买的服务器还在路上。就在这走投无路的情况下,我突然想起2个月之前和这个设备的负责人说要去测试一下(换句话说,我放了他两个月的鸽子。。) 看完软件说明书,我就坐着地铁揣着硬盘,硬盘里装着一个260M的基因组和230个GBS测序的数据(80G)跑到仁科生物公司以测试软件之名实为蹭别人的服务器。 首先我把数据一股脑地全从硬盘里拷到固态硬盘挂载的 /staging下 然后是建立索引: dragen -h-ht-reference reference.fa --output-directory reference --build-hash-table true 接着我现场写了一个shell脚本用来批量分析,命名为 run_dragen.sh #!/bin/bash set -e set -u set -o pipefail REF=$1 SAMPLES=$2 samples=$(cat $SAMPLES) # loading hash table dragen -l -r $REF # calling gvcf for each sample mkdir -p GVCF for sample in ${samples} prefix=$(basename ${sample}) if [ ! -f GVCF/${prefix}.done ]; then dragen -f -r $REF \ -1 ${sample}.1.fq.gz -2 ${sample}.2.fq.gz \ --enable-variant-caller true \ --vc-emit-ref-confidence GVCF \ --vc-sample-name ${prefix} \ --output-directory GVCF \ --output-file-prefix ${prefix} \ --enable-duplicate-marking false \ --enable-map-align-output true touch GVCF/${prefix}.done find GVCF/ -name "*.gvcf.gz" | grep -v "hard" > gvcfs.list # merge gvcf and join calling mkdir -p vcf_result if [ -f gvcfs.list -a ! -f vcf_result/combine.gvcf.gz]; then dragen -f -r $REF \ --enable-combinegvcfs true \ --output-directory vcf_result \ --output-file-prefix combine\ --variant-list gvcfs.listfi if [ ! -f vcf_result/join_calling.vcf.gz ] dragen -f -r $REF \ --enable-joint-genotyping true \ --output-directory vcf_result \ --output-file-prefix joint_calling\ --variant vcf_result/combine.gvcf.gzfi 运行方法: # 创建一个文件存放待分析的样本 find /staging/xuzhougeng/00-raw-data/ -name "*.fq.gz"| sed 's/\(.*\)\.[12].fq.gz/\1/' | uniq > samples.txt # 执行命令, 参数分别是索引的文件夹和样本文件 bash run_dragen.sh reference samples.txt &> run.log & 按照我的估算,每个样本至少得要花个20分钟得到GVCF文件吧,毕竟我用BWA-MEM10个线程进行比对也要10min呀。事实证明我还是低估了程序猿的能力值,每个GBS样品得到GVCF文件居然只要不到1min。。 RUN TIME,,Total runtime,00:00:56.528,56.53 得到的GVCF可以进行合并,但是有一个问题,就是超过200样本就会出错,而且Join calling运行也不需要combine,所以后续的代码就删掉了merge这一步 ...find GVCF/ -name "*.gvcf.gz" | grep -v "hard" > gvcfs.list # join callingif [ ! -f vcf_result/join_calling.vcf.gz ]t hen dragen -f -r $REF \ --enable-joint-genotyping true \ --output-directory vcf_result \ --output-file-prefix joint_calling\ --variant-list gvcfs.listfi 在Joint Calling步骤花的时间比较长,时间是"RUN TIME,,Total runtime,00:14:24.272,864.27". 虽然软件运行速度是很快,但是写出上面的代码并且调试却花了我好久时间,于是这两天时间我就在公司里敲代码。除了GBS数据,第二天我还带着另一个260M基因组(Canu初步组装和arrow polish后到版本)和一个100X重测序数据(压缩后10G数据)去测试,分别在固态硬盘和我的移动硬盘里测试,结果如下: 固态硬盘 draft.fa: 06:31(mapping) + 28.17 (varaint calling) 普通移动硬盘 draft.fa: 28.33(mapping) + 29.03 (varaint calling) 固态硬盘 polish.fa: 06:09(mapping) + 13.093(variant calling) 普通移动硬盘 polish.fa: 27.46(mapping) + 13.33(varint calling) 这里有两个结论 是比对的IO对速度影响非常大,也就是要一定要在固态硬盘里发挥它最大的威力。 重测序样本与参考基因组的差异程序影响variant calling这一步。 从上面的测试而言,DRAGEN的运算速度的确是非常快的。虽然你需要先把拷贝数据这一步会花点时间,但是你从公司拿到的数据其实也要拷贝到服务器才行,所以拷贝数据是不可避免的。 对于公司而言,原本需要两天才能跑完的分析可能现在2小时或者1小时不到就能搞定了,那么业务速度就快了,此外也不需要搭建spark或者自己搞一套对GATK进行并行,更何况GATK商用是要钱的,国内很多公司都是偷偷的在用吧。对于医院而言,嗯,他们不差钱。对于科研机构而言,除非专门搞一个平台管理,不然一年花掉2T的分析数据量还是有难度。 PS:当然都是比对,所以这个软件也能用于分析RNA-seq,ChIP-seq,ATAC-seq等illumina高通量测序数据,但是三代测序数据目前搞不定,不知道未来会不会支持。 PPS: EDICO公司看到这篇文章后请给我打广告费

LAI: 评估基因组质量一个标准

基因组组装完成之后,就需要对最后的质量进行评估。我们希望得到的contig文件中,每个contig都能足够的长,能够有一个完整的基因结构,归纳一下就是3C原则: 连续性(Contiguity): 得到的contig要足够的长 正确性(Correctness): 组装的contig错误率要低 完整性(Completeness):尽可能包含整个原始序列 但是这三条原则其实是相互矛盾的,连续性越高,就意味着要处理更多的模糊节点,会导致整体错误率上升,为了保证完全的正确,那么就会导致contig非常的零碎。此外,这三条原则也比较定性,我们需要更加定量的数值衡量,目前比较常用的标准是N50和BUSCO/CEGMA。 最近有一篇文章"Assessing genome assembly quality using the LTR Assembly Index (LAI) "提出用长末端重复序列来评估基因组完整度,因为LTR比较难以组装,于是就用作评估结果的一个参数了。那问题来了,什么是LTR序列,LTR是在原病毒(整合的反转录病毒)两末的重复序列,结构见下图 LTR结构 上图中TSD表示target site duplications,红色三角表示LTR motif。A图是一个完整的LTR结构,其中a,b,c是LTR_retriever的分析目标。 LAI指数就是完整LTR反转座子序列占总LTR序列长度的比值。 其实作为一个农学出身,看到LAI,我脑海就想到了Leaf Area Index(叶面积指数) 本文以拟南芥的基因组为例来测试一下这个软件 要想保证软件能够顺利的安装,需要先安装如下这几个软件, 好消息是这些软件都可以通过bioconda解决 makeblastdb, blastn, blastx cd-hit-est hmmserch RepeatMasker 然后从GitHub上下载软件 cd ~/opt/biosoft git clone https://github.com/oushujun/LTR_retriever.git 进入LTR_retriever文件下修改paths文件,提供每个软件所在的文件路径,下面是我的配置,你需要按照实际所在路径来设置 BLAST+=/home/xuzhougeng/opt/biosoft/ncbi-blast-2.7.1+/bin/ RepeatMasker=/home/xuzhougeng/opt/biosoft/RepeatMasker/ HMMER=/home/xuzhougeng/opt/anaconda2/envs/maker/bin/ CDHIT=/home/xuzhougeng/opt/anaconda2/envs/assembly/bin/ 此外,你还需要安装GenomeTools或者LTR_FINDER,或者MGEScan_LTR才能提取出LTR序列,我这里下载的是LTR_FINDER cd ~/opt/biosoft git clone https://github.com/xzhub/LTR_Finder.git cd LTR_Finder/source/ 第一步让我们用LTR_FINDER找到基因组的LTR序列 ~/opt/biosoft/LTR_Finder/source/ltr_finder -D 20000 -d 1000 -L 700 -l 100 -p 20 -C -M 0.9 Athaliana.fa >Athaliana.finder.scn 这里的-D表示5'和3'LTR之间的最大距离,-d表示5'和3'LTR之间的最小距离,-L表示5'和3'LTR序列的最大长度,-l表示5'和3'LTR序列的最小长度,-p表示完全匹配配对的最小长度,-C表示检测中心粒(centriole)删除高度重复区域,-M表示最小的LTR相似度。如果不怎么该怎么设置就用默认值。 第二步运行LTR_retriever根据LTR_FINDER的输出识别LTR-RT,生成非冗余LTR-RT文库,可用于基因组注释 ~/opt/biosoft/LTR_retriever/LTR_retriever -threads 4 -genome Athaliana.fa -infinder Athaliana.finder.scn 这里的-infinder表示输入来自于LTR_FINDER,它支持同时输入LTRharvest的输出(-inharvest)和 MGEScan-LTR 的输出(-inmgescan). 嫌速度太慢,可以用-threads增加线程数 这一步会调用RepeatMasker,而RepeatMasker要求序列ID长度不大于50个字符,所以请在第一步的时候请先对ID进行修改。 第三步计算LAI。如果前面找到LTR序列太少,低于5%,这一步程序就会报错,那么你就需要调整第一步参数,可能是太严格了。 /opt/biosoft/LTR_retriever/LAI -t 10 -genome Athaliana.fa -intact Athaliana.fa.pass.list -all Athaliana.fa.out 这里最后的结果文件为Athaliana.fa.out.LAI, 第二行就是总体信息,其中RAW_LAI是12.88, LAI是14.47 Chr From To Intact Total raw_LAI LAI whole_genome 1 119667750 0.0079 0.0612 12.88 14.47 得到的LAI值按照如下评估标准进行分类: Category Examples Draft 0 ≤ LAI < 10 Apple (v1.0), Cacao (v1.0) Reference 10 ≤ LAI < 20 Arabidopsis (TAIR10), Grape (12X) 20 ≤ LAI Rice (MSUv7), Maize (B73 v4) 和例子一样,TAIR10是中等水平。 参考文献: Ou S. and Jiang N. (2018). LTR_retriever: A Highly Accurate and Sensitive Program for Identification of Long Terminal Repeat Retrotransposons. Plant Physiol. 176(2): 1410-1422. Ou S., Chen J. and Jiang N. (2018). Assessing genome assembly quality using the LTR Assembly Index (LAI). Nucleic Acids Res. gky730: https://doi.org/10.1093/nar/gky730

学习从来都不是一件容易的事

不知不觉间,我给别人上了好多次课程,有些是线上课程,有些是线下课程,有些是收费的,有些是免费的,如果有什么感悟要说,那就是“学习从来都不是一件容易的事。 思考一个问题,你是如何学会数学的乘法和除法的?为什么我要为这个问题呢,这要从这次回家休假,我接到了一个非常高难度的任务说起。我给自己的侄女辅导数学,其中有一个题目,“现在有三个牌堆,分别有6,7,8张扑克牌,如何让这三个牌堆的有一样多的牌呢?” 你可能会说很简单呀,一共21张,平均每堆7张,所以把8张移动1张给6张就行了呀。我相信你如果给才学会用手指头进行加减运算的孩子这样讲,他肯定不知道你在说什么,一脸问号的看着你。然后你就开始不耐烦,说,怎么那么简单的运算你怎么就不能理解呢? 那些你理所当然的事情,其实你花了大量的时间进行练习才掌握。 但其实当年的你也是那么痛苦,学会你现在认为“非常简单”的运算,你忘记了那个状态,想不明白为什么要那么复杂。 最后,我能想到的一个方案就是,给她找了21张纸,分成6,7,8三堆,让她重新分配,分成7,7,7,然后问她前后发生了什么变化,问她把21平均分配3个人是多少张,最后才和她说这就是除法。当然即便和她说了“除法”这个概念,她依旧不知道“除法“到底是什么。她还要不断做题目,总结规律,或许才能学会了除法。或许再过几个月,大脑又发育了一些,就突然学会了。 因此,学习本身就非常的困难,每掌握一个技能都要耗费一定量的精力,然后在这个基础上才能学习更高级一些概念,如此这般才能循序渐进,渐入佳境。 但是在这个贩卖焦虑的时代,很多人似乎忘记了自己掌握目前技能所耗费的精力,以为自己只要花一些钱,就能身心愉悦地把技能掌握了。花了几千块钱上线下课,全程不想费脑子,一旦遇到问题就寻求他人的帮助,而不是自己尝试解决。最后自己一无所获,除了几个能用来给别人吹牛的概念。 可是学习过程从来不是愉悦的,毕竟开心的瞬间转眼即逝,痛苦的记忆在你脑海里才会深刻。 我目前学习就是这个状态,每天都得去接触陌生的领域,硬着头皮探索,做笔记记录自己目前的心得,然后通过实践巩固,大部分的技能都是这样子慢慢学会。很多人说我学习速度快,其实只是领域不同,我必须要花足够多的时间才能让自己混口饭吃。

三年后,我又来到了北京

我再次来到了北京,距离上次,时间已过去3年 三年前的那个暑假,机缘巧合,我去参加北京植物所开的光合夏令营。当时来回只能报销火车的硬座票,所以我硬生生从发车的下午坐到了隔天的上午,坐着实在是不舒服,一宿未眠。去的14小时的时间我看完了几百页的爱因斯坦传,回来的14小时我看完了富兰克林传。当然,当年看书主要是打发时间,所以早已忘记那些书到底讲了什么,也不知道对我的人生轨迹有多少影响。 在北京夏令营的时间里,住的是香山会议中心,可比另一批住如家的人高级多了。那时候认识了一些小伙伴,跟他们坐公交车去了北京的几个地方。南锣鼓巷里全是人和吃的,鸟巢和水立方夜景很美,中国农大去看了一眼实验室,清华大学根本进不去只好悻悻而归。 三年前的那些人,那些事就这样存在我的记忆里和照片里。那时的我也无法知晓我会在三年后会再次来到北京,还会去见那些熟悉而陌生的人。 这次去北京也有点坎坷,台风将我困在了上海一天,一度要放弃计划,也打乱了我原先的所有的计划,但最终还是去了。课程就让我略去不谈吧,尽管此行的主要任务是培训,让我写点流水账说说那些人。 在一只烤鸭的诱惑下,我又去中国农业大学逛了一小圈,也终于在天未黑时在农大老校门拍照打卡,可惜的是烤鸭只能吃半只,因为肚子能承受的有限。随后又在另一个小伙伴的掩护下顺利进入了清华大学,只可惜天色已黑加上阵雨袭来,匆匆逛了校园就回去了。如有可能,真希望能在这里学习生活一段时间,而不是作为游人匆匆来过。最后又回到了中农大看了国内仅有的一批上个世纪九十年代繁衍到现在的虫子,它们见证了一个实验室的几代院士。 另一晚去见了当年的室友,本想着去他的学校一边逛逛再随便聊聊,但最后见面的地方是一家火锅店,并且还带着他的妻子秀了我一脸。当年他考上北林之后,他的未婚妻也决定来北京找工作和他一起继续生活,所以在去年终于领了证,是我们寝室最早成家的人。不知为何,我想起了大学时在实验室结识的一个人,再次翻看他的朋友圈发现封面已是婚纱照,而几句交流之后了解到半年后他就可以晒娃了。 结束培训的次日,北京又恢复了它的燥热。 2017年的11月份我曾写过一片文章关于果壳生物公司的招聘,那时候办公桌空荡荡,虚位以待,而我在参观时发现他们团队虽已有了20多人,但还是满负荷运算,听他们说还需要继续扩大团队。也不知道在我毕业那年,他们公司会发展的如何呢? 匆匆参观完果壳生物,就赶去中科院的遗传发育所去见另一位小伙伴,同时还幸运的见到了公众号《宏基因组》的创始人刘永鑫。第一眼,丝毫无法把眼前的他和微信号头像联系起来,他也不知道我是谁,直到我说出“生信媛”。他带我参观他们所曾组装了小麦基因组的服务器,想象下小麦的基因组的组装难度,就知道这台服务器的算力有多强大,而我自己所里的服务器平台却无人管理。原本打算稍微看看他们的实验室就可以离开了,但是他们所里的顶层居然还有咖啡厅,于是我们就一边喝着咖啡,一边谈饱和的公众号、听他列举生信届的大牛。 北京的最后一站是陈同老师的易汉博公司。和陈老师聊的更久,吐槽生信培训中的那些坎坷,听他讲述在科研服务模式上的探索。很多事情说出来的时候轻描淡写,但背后的努力其实超乎想象,愿他的公司能够稳步发展吧。 脑中有很多话想说,而写下来的却只这些,情深词乏,不知如何表达。

一文了解Perl语言

我在公众号发过很多编程语言的学习笔记,但是一直没有发Perl语言的编程教程。我大学的时候,学过一段时间的Perl语言,所以和Perl也有点缘分。这次去北京参加培训时发现他们教的Perl,所以接着机会发一波我现场的学习记录。 什么是Perl 学习一门语言最好了解下它的历史,知道它能干什么,有什么优势是,有什么不足。以下内容来自于百度百科 Perl,一种功能丰富的计算机程序语言,运行在超过100种计算机平台上,适用广泛,从大型机到便携设备,从快速原型创建到大规模可扩展开发。 [1] Perl最初的设计者为拉里·沃尔(Larry Wall),于1987年12月18日发表。现在的版本为Perl 6,于2015年12月25日更新。 Perl借取了C、sed、awk、shell 脚本语言以及很多其他程序语言的特性,其中最重要的特性是它内部集成了正则表达式的功能,以及巨大的第三方代码库CPAN。简而言之,Perl像C一样强大,像awk、sed等脚本描述语言一样方便,被Perl语言爱好者称之为“一种拥有各种语言功能的梦幻脚本语言”、“Unix 中的王牌工具”。 Perl 一般被称为“实用报表提取语言”(Practical Extraction and Report Language),你也可能看到“perl”,所有的字母都是小写的。一般,“Perl”,有大写的 P,是指语言本身,而“perl”,小写的 p,是指程序运行的解释器。 从这段描述中,我们就知道了Perl最重要的就是它的文本处理能力, 而早期生物信息学的本质就是序列分析,所以“确认过眼神,我遇见对的人”,Perl成了生物信息必学语言。 Perl的安装 Perl语言在Linux系统以及其他类Unix系统,比如说MacOS里都是内置的,所以不需要额外安装,而在Windows系统中则需要额外下载,下载站点当然是官方的<www.perl.org> 我后续在MacOS上学习Perl, 版本是5.18.2,这版本有点老旧,但是不妨碍学习。 Perl的Hello World 大部分的教程都是让大家写一个Perl脚本,比如说hello.pl,代码如下 #!/usr/bin/env perl print "hello world \n" 然后用perl hello.perl执行。当然最快的方法还是用Perl的一行命令, 如下所示。 perl -e 'print "hello world"' 据说Perl高手能把Perl一行命令用的出神入化 Perl的编程基础 你可以用20多门语言说“我爱你”,但是不代表你会20门语言,所以会用Perl写"hello world"其实只是开始而已。了解一门编程语言,最少要掌握如下内容: 数据类型和数据结构 控制语句:条件语句、循环语句的用法 函数写法和常用函数 如何调用其他包 数据类型和数据结构 Perl支持的数据类型其实和大多数编程语言差不多,都有数值型(整数和浮点数)、字符串 数据结构分为:标量、数组、哈希。 这些数据结构都会和变量联系在一起,才能根据变量调取数据。Perl和Python一大不同就在于,Python认为括号",'以外的字符串都是变量名,而Perl的变量名之前必须要有专门的符号指明,这和shell类似。 perl -e '$id = abc; print "hello $id"' ; # 标量定义和调用 perl -e '@number = (1,2,3,4); print "@number\n"'# 数组定义和调用 perl -e '@number = (1,2,3,4); print "$number[1]\n"'# 数组定义和调用 perl -e '$age{Tom}=25;$age{David}=10; print "$age{Tom}\n"' # 哈希定义和调用 perl -e '%hash = (a=>b,c=>d,e=>f); print "$hash{a}\n"' # 哈希定义和调用 所以以后在Perl脚本看到一堆的符号时不要慌,看到$是标量,看到@是数组, 看到%就是哈希。 哈希对应Python的字典(dict), 对应R的列表(list) 据说所有逻辑都可以用条件语句和循环语句来编写,所以每一门编程语言都一定要有。 条件语句的写法 perl -e '$x=0; if($x>0){$y=10;}elsif($x<0){$y=-10;}else{$y=0;} print "$y\n";' 循环语句有三种写法,for ,foreach, while for ($x=5;$x<10;$x++){ print "$x\n"; foreach $x(1..5){ print "$x\n"; while 比较适用于文件读取和输出, 让我们先创建一个文件cat /etc/passwd > passwd.txt, 然后输出到另一个文件中 #/usr/bin/env perl my $file_in = "passwd.txt"; my $file_out = "passwd_out.txt"; open IN, "<", $file_in; open OUT, ">", $file_out; while (<IN>){ print OUT "$_\n"; close IN; close OUT; 这里出现几个陌生的东西,一个是my表示声明变量,取决于具体位置来确定定义的变量是局部变量还是全局变量。 <>表示读取文件 , open是一个函数,用于创建一个文件句柄从而进行文件读写, close就是关闭文件读写。 这里有一个诡异的符号是的$_, 它表示是读取的当前行 函数写法和常用函数 在Perl里是找不到对象的,因为它不是面向对象编程的语言,它只有函数。 # 声明函数 sub sum{ my ($m, $n) = @_; return ($m + $n ); # 调用函数 $sum_number = &sum(2,8); print "$sum_number \n" 这里又出现了一个诡异的字符@_ ,表示是函数传参时放置参数的数组,可以从中取实参。 调用函数的时候需要在函数名前用&声明。 当然,大部分情况下,我们是不怎么写函数,我们都是面向调用函数编程,所以知道Perl里有哪些好用的函数很重要 字符串相关函数: split, join, substr, chomp 数组相关函数:shift, pop, unshift, push, scalar, length 哈希相关函数:keys,values, exists -其他函数: sort, reverse 关于函数的用法,在http://perldoc.perl.org/perl.html进行搜索,一般要了解一门编程语言一定要查不知道多少遍的帮助文档。 正则表达式 其实正则表达式才是Perl的重点,学Perl不用Perl的正则表达式等于是没有学过Perl,但是这部分内容其实是非常多的,这里举出一个例子, $seq = "ACCGGATCATTGTCAA"; if ($seq =~ m/AC/){ print "$seq \n"; $seq =~ s/AA/GG/; print "$seq \n"; $seq =~ tr/ATGC/TACG/; print "$seq \n"; 看到的=~就意味着后面将会出现正则表达式。=~不是简单的赋值符号=,也不是==用于逻辑判断,而是将变量的值进行对应的正则运算 m/AC/ : 指的是匹配(match),和命令行的grep一样 s/AA/GG/: 指的是替换(substitute), 将符合要求的替换成对应的值 tr/ATGC/TACG/: 指的是转换,和命令行的tr一样 希望本文让你对Perl有一个大致印象,想要继续深入的话,推荐两本书: Perl语言入门:俗称小骆驼 精通正则表达式:真正讲透了正则的精髓