文献笔记三十三：结合二代三代测序数据组装叶绿体基因组

文章

Assembly of chloroplast genomes with long- and short-read data: a comparison of approaches using Eucalyptus pauciflora as a test case
2018 BMC Genomics
Australian National University

研究内容

以Eucalyptus pauciflora为例，探索组装叶绿体基因组最有效的方法

DNA提取测序

二代测序
leaves
total DNA
150bp paired-end sequencing with a roughly 400 bp insert size
质控
三代测序
High molecular weight DNA
Libraries were prepared according to the ONT 1D ligation library protocol.
质控

组装方法

二代测序数据
unicycler软件
三代测序数据
Hinge软件： an assembler designed for solving the long repeats problem in long-read assemblies of circular genomes
Canu软件
结合二代三代数据
unicycler软件

软件安装

unicycler

## 创建虚拟环境
conda create -n chloroAssembly python=3.6
conda activate chloroAssembly
conda install unicycler
###删除虚拟环境 conda remove -n chloroAssembly --all

试着运行unicycler软件主页中的例子

https://github.com/rrwick/Unicycler
使用的数据可以在软件主页找到下载链接

unicycler -1 short_reads_1.fastq -2 short_reads_2.fastq -l long_reads_high_depth.fastq -o output_dir -t 16

数据是Helicobacter pylori，在NCBI查了一下基因组大小1,667,867bp，使用unicycler的组装结果

grep ">" assembly.fasta
>1 length=1645796 depth=1.00x circular=true

稍微有点差别，可能是不同的株系吧我猜

graph.png

下载Eucalyptus pauciflora的全基因组测序数据

三代测序数据

wget ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR715/SRR7153095/SRR7153095.sra
fasterq-dump SRR7153095.sra -p

检查三代测序数据质量，使用到的是fastqc软件，原来fastqc软件还可以用于三代测序数据

mkdir qcResult
fastqc SRR7153095.sra.fastq -o qcResult -t 8

去除接头，用到的软件是porechop
https://github.com/rrwick/Porechop

conda install porechop
porechop -i SRR7153095.sra.fastq -o longReadsRemoveAdapter.fastq -t 8

数据过滤，质量值大于9，最小长度5000，使用到的软件是nanofilt

conda install nanofilt
bgzip longReadsRemoveAdapter.fastq
zcat longReadsRemoveAdapter.fastq.gz | NanoFilt -q 9 -l 5000 > longReadsRemoveAdapterTrim.fastq

二代测序数据

数据下载

wget ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR715/SRR7153063/SRR7153063.sra
fasterq-dump --split-files SRR7153063.sra -p
wget ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR715/SRR7153071/SRR7153071.sra
fasterq-dump --split-files SRR7153071.sra -p

数据过滤，不按照论文中提供的脚本来了，直接使用fastq软件进行过滤了
软件主页 https://github.com/OpenGene/fastp

fastp -i SRR7153071.sra_1.fastq -I SRR7153071.sra_2.fastq -o shortReads71_R1.fastq -O shortReads71_R2.fastq
fastp -i SRR7153063.sra_1.fastq -I SRR7153063.sra_2.fastq -o shortReads63_R1.fastq -O shortReads63_R2.fastq

未完待续 ......

20200105继续

数据预处理做完了，接下来是从总DNA里提取数据叶绿体基因组的reads

论文里使用到31个Eucalyptus叶绿体基因组作为参考序列，我这里直接使用论文的组装结果一条序列作为参考，另外论文里还写到

The chloroplast genome is circular, but alignment algorithms rely on linear genomes. Simple linearization of the reference set would risk failing to capture reads that span the point at which the genomes were circularized. To avoid this, we duplicated and concatenated the sequence of each genome in the reference set.

论文中的组装结果也可以直接在论文里提供的github主页上下载到

三代测序数据筛选叶绿体基因组reads，使用到的软件是 blasr

conda install blasr
blasr longReadsRemoveAdapterTrim.fastq ../reference/epauCPgenome.fasta --nproc 8 --bestn 1 -m 1 --minMatch 15 --minAlnLength 5000 --out longReads.blasr.out
conda install seqtk
awk '{print $1}' longReads.blasr.out > subsetid.txt
seqtk subseq longReadsRemoveAdapterTrim.fastq subsetid.txt > longReadsCPgenome.fastq

二代测序数据提取叶绿体基因组reads，使用到的软件是bowtie2

bowtie2-build epauCPgenome.fasta epau
bowtie2 -q -x epau -1 ../short_reads/shortReads63_R1.fastq -2 ../short_reads/shortReads63_R2.fastq -p 8 -S shortReads63.sam
samtools view -S -b -o shortReads63.bam shortReads63.sam
samtools sort -@ 4 -O bam -o shortReads63.sorted.bam shortReads63.bam
samtools index shortReads63.sorted.bam
samtools view -u -f 1 -F 12 shortReads63.sorted.bam > shortReads63.sorted.aligned.bam
conda install bedtools
samtools sort -n -O bam -o shortReads63.sortedbyname.aligned.bam shortReads63.sorted.aligned.bam ###-n参数sort bam文件按照reads的名字
bamToFastq -i shortReads63.sortedbyname.aligned.bam -fq mappedcpR1.fastq -fq2 mappedcpR2.fastq

数据组装

单纯用三代测序数据

canu -p dd -d onlylongreads genomesize=160k -nanopore-raw longReadsCPgenome.fastq

组合二代三代测序数据

unicycler -1 ../reference/mappedcpR1.fastq -2 ../reference/mappedcpR2.fastq -l ../long_reads/longReadsCPgenome.fastq -o output_dir -t 16

结果

grep ">" assembly.fasta
>1 length=159945 depth=1.00x circular=true

论文中提到的最终组装长度为159942bp.
欢迎大家关注我的公众号
小明的数据分析笔记本

公众号二维码.jpg