Pre-release! 在上一篇中准备好了需要的二代测序数据，这篇就开始call SNP bwa比对 构建索引 在开始比对之前要先对参考基因组构建索引。需要使用到三个软件sentioen bwa（目前看来跟bwa本身没有太多区别），samtools和picard，这里主要介绍samtools和picar的安装。 # samtools conda install -c bioconda samtools # picard wget https://github.com/broadinstitute/picard/releases/download/2.25.6/picard.jar java -jar picar.jar 构建索引的流程也很简单bwa index生成.amb，.ann，.bwt，.pac和.sa文件，samtools faidx生成.fai，picard CreateSequenceDictionary生成.dict sentieon bwa index ${GENOME} samtools faidx ${GENOME} java -jar picard.jar CreateSequenceDictionary REFERENCE=${GENOME} OUTPUT=${DICT} 比对与排序 这里将比对和排序两个步骤合并在一起执行。比对的部分还是使用的sentioen bwa，算法部分使用的是mem，-M ( sentieon bwa mem -M -R "@RG\tID:${RGID}\tSM:${SM}\tPL:${PL}" -t ${NT} -K 10000000 ${GENOME} ${FQ1} ${FQ2} || echo -n 'error' ) | sentieon util sort -r ${GENOME} -o ${SORTED.BAM} -t ${NT} --sam2bam -i - 统计信息 去除重复 Haplotyper

背景 群体分析通过对多个个体进行重测序，可以发现不同个体中存在的变异，挖掘在驯化过程中收到选择的位点。目前网上大多的教程都是基于GATK进行的call SNP，但是GATK对多线程的支持很不友好，在GATK3中可以通过-nt和-ncr参数来进行多线程运算，在GATK4中已经取消了这些参数，虽然还是可以通过调整–native-pair-hmm-threads来加快核心计算程序，但是效果很不理想。相较之下Sentieon能提供更稳定更快的计算，因为是基于GATK的，结果也较为准确。Sentieon提供了很详细的帮助文档，我在此基础上记录下操作过程可供参考。 这个系列大致分为三部分吧，第一部分是数据的准备，第二部分是call SNP，第三部分是其他下游分析。 FastQC质控 测序公司在交付重测序原始数据的时候基本已经通过了质量控制，交付的是clean data，而且现在二代测序十分普及，流程很完善，基本可以直接使用公司提供的原始数据了。这一步主要是针对在NCBI上存储的其他研究测的数据，特别是早些年的数据质量参差不齐，需要自己质控一下。这里使用的软件是FastQC，可以通过conda直接安装。 conda install -c bioconda fastqc 如果想一次对多个文件提交操作，可以使用下面的命令，此处-t设置线程数为108是因为有108个fq文件，每个文件使用一个线程并行计算。-o指定输出的文件路径，此文件夹需提前建立。 fastqc -t 108 -o fastqc/ *.gz FastQC会为每个fq文件生成一份.html报告和一个.zip压缩包，在本地电脑浏览器打开.html文件就可以查看报告。但是如果文件数量太多一个一个看太不方便，可以使用MultiQC软件汇总结果，生成一份总的报告。MultiQC可以用过pip安装。 pip3 install multiqc 直接提供fastqc的结果文件夹给MultiQC即可，它会检索文件夹里报告的数量，生成一份汇总报告。 multiqc fastqc/ 结果文件是一个.multiqc_report.html和一个multiqc_data文件夹，下载html文件在浏览器就可以打开查看报告了 Trimmomatic过滤 对于原始数据质量不过关的文件，可以使用Trimmomatic进行过滤。主要的过滤目的有两个，一个是去除低质量的碱基，还有一个就是去除接头。Trimmomatic使用的是java，在官网下载解压就可以用。 wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.39.zip unzip Trimmomatic-0.39.zip cd Trimmomatic-0.39 java -jar trimmomatic-0.39.jar 这里介绍使用的几个参数。PE是双端测序模式，如果是单端测序改成SE。-threads指定线程数量，我试过给一个比较大的数，但是cpu的使用率一直在300多一些，所以每个文件就指定了4个线程。后面的文件名顺序依次是原始fq1 原始fq2 过滤后fq1 过滤后未成对fq1 过滤后fq2 过滤后未成对fq2。ILLUMINACLIP指定接头文件的路径和匹配参数，这里的匹配参数具体意义还不明确，暂时使用陈连福老师教程里的数值。LEADING指定去除前端碱基质量值(Qscore)小于20的碱基。TRAILING指定去除末端碱基质量值(Qscore)小于20的碱基。SLIDINGWINDOW指定以4个碱基为一个窗口滑动，去除平均质量值小于25的碱基。MINLEN指定保留的最小读段长度。TOPHRED33将碱基质量值转换为phred33格式（原始就是phred33格式的话应该可以不加，但是暂时还是按照陈连福老师的教程加了）。 java -jar ~/software/Trimmomatic-0.39/trimmomatic-0.39.jar PE -threads 4 ${RAW.1.FQ.GZ} ${RAW.2.FQ.GZ} ${CLEAN.1.FQ.GZ} ${UNPAIR.1.FQ.GZ} ${CLEAN.2.FQ.GZ} ${UNPAIR.2.FQ.GZ} ILLUMINACLIP:${ADEPTER.FA}:2:30:10 LEADING:20 TRAILING:20 SLIDINGWINDOW:4:25 MINLEN:50 TOPHRED33 对于多个文件都要过滤的话可以参考下面的脚本模板一次性提交...

背景 前段时间开始了解并使用GitHub，看到大佬们自己在博客上分享帖子，很是羡慕，借此机会打算利用GitHub Pages来搭建一个个人网址。搭建的软件（暂且这么称呼吧）比较常见的有GitHub自己推荐的基于Ruby的Jekyll，Heko和基于GO的Hugo。简单尝试了Jekyll之后被复杂的配置方式劝退了，恰巧看到v站上关于Heko和Hugo的讨论，于是决定试一下Hugo（但还是要夸一下Heko优秀的中文文档），以下主要是记录流程，大部分都是参考网上现有的教程，但是找不到参考教程的网址了，实在抱歉，后面找到会再补充到文末。 这个系列文章应该会分成三篇，第一篇介绍Hugo的安装与使用，第二篇介绍GitHub Pages的申请与配置，第三篇讲利用PicGo和sm.ms实现图床功能。 Hugo的安装与使用 这部分主要参考的是Hugo中文文档，虽然好像已经很久没有更新了，但是还是顺利安装了下来。因为我使用的是macOS的操作系统，所以直接通过brew安装即可。 brew install hugo 安装完成之后就可以通过hugo --help查看Hugo的详细参数啦。接下来开始创建一个站点，站点保存的路径我用PATH替代，这里可以根据需要换成自己的路径，可以是相对路径也可以是绝对路径。（类似的方式在下文中会经常出现，就不再一一交代，即大写字母替换为自定义的参数） hugo new site PATH 这样在你的路径下就生成了一个文件夹，文件夹内保存着这个站点的所有文件。其中，content文件夹主要存放后续的文章，theme文件夹存放主题文件，config.tmol是设置站点的样式。Hugo提供了非常多的预设主题，可以直接在主题商店下载使用。这里我选了一个叫PaperMod的主题，复制下载链接，在theme文件夹内通过git clone下载即可。 git clone https://github.com/adityatelange/hugo-PaperMod.git 之后回到config.tmol中添加一句theme = "hugo-PaperMod"就可以应用这个主题了。顺带提一句，这个文件中的title就是站点的名称，可以更改为自己想要的名字。 设置好了站点主题，就开始新建一篇文章吧，默认直接采用markdown的后缀名。 hugo new post/POST_NAME.md 这样在content/post文件夹下就会生成一个.md文件，打开里面会包含由---分隔的头信息，其中title就是文章的标题，draft表示文章的状态，默认一开始就是草稿的状态，这个时候在hugo server -D预览的时候能顺利显示该文章，但是在hugo发布文章的时候并不会显示，所以当你要发布的时候记得把draft改为false。 写好了文章就在本地浏览器里预览效果吧。回到站点文件夹，运行hugo server -D之后，在浏览器内访问地址http://localhost:1313就可以看到显示的效果了。如果一切满意，就可以直接运行hugo（对，没有错，不用加任何参数，直接hugo）发布站点。这时候会生成一个public的文件夹，这个文件夹内就包含了由hugo生成的html文件，之后通过把这些文件同步到GitHub Pages上就可以实现一个个人网址啦。

Pre-release! Prepare input data 准备输入文件 GFF cds Format 对输入文件进行统一的格式化 convert GFF to BED 将GFF格式转换成BED格式 python3 -m jcvi.formats.gff bed --type=mRNA --key=Parent JamaicanLionDASH.gene.gff3.gz -o CsJLD.bed formatting cds sequences python3 -m jcvi.formats.fasta format JamaicanLionDASH.cds.fasta.gz CsJLD.cds Pairwise synteny search 搜索共线性对 python3 -m jcvi.compara.catalog ortholog CsFN CsPK --no_strip_names 此步骤会调用LAST以CsPK为subject建库,以CsFN为query做blastp。结果文件中.last是比对的全部结果，.last.filtered是过滤后的结果，但过滤的规则尚不清楚。同时生成的.anchors和.lifted.anchors文件格式尚不清楚。生成的.pdf文件绘制了CsFN和CsPK的共线性点图。 若杂点太多可以通过添加--cscore=.99参数修改last的C-score过滤阈值，在此之前需删除旧的.last.filtered文件。 Macrosynteny visualization 宏观共线性可视化 seqids and layout chr1,chr2,chr3,chr4,chr5,chr6,chr7,chr8 Pp01,Pp02,Pp03,Pp04,Pp05,Pp06,Pp07,Pp08 seqids文件是CsFN和CsPK中需要展示的染色体编号，编号的排列顺序与结果文件中染色体的排列顺序一致 两个个体之间的染色体编号若有重复，可能导致在微观共线性可视化时产生错误。 # y, xstart, xend, rotation, color, label, va, bed ....