Plant contig clustering based assembly (PLCL) pipeline: an efficient long-read assembly tool for plant chloroplast and mitochondrial genomes

sratoolkit.2.11.0-ubuntu64/bin/fastq-dump ERR2173373 –gzip

NanoPlot \
	-t 8 \
	--fastq $HOME/PLCL_At / ERR2173373.fastq.gz \
	-o OUTPUT_dir \
	--loglength

NanoPlot \
        -t 8 \
        --fastq $HOME/ONT_genome_raw.fastq.gz \
        -o OUTPUT_dir_raw \
        --raw \
        --store

library(ggplot2)
df <- read.table("NanoPlot-data.tsv", header = T)
ONT <- ggplot(df, aes(x=lengths)) + geom_histogram(bins=1000) +
  scale_x_log10(limits=c(1,1000000),labels=scales::comma) +
  ylim(0,10000) +
  xlab("Read_length") +
  ylab("Number_of_reads")+
  theme(text=element_text(size=20))+
  geom_vline(xintercept = c(3500), linetype="dashed", color="yellow")
ggsave("ONT.emf",plot=ONT)

wtdbg2 \
        -x ont \
        -g 135m \
        -i $HOME/PLCL_At / ERR2173373.fastq.gz \
        -t 32 \
        -fo Atont

wtpoa-cns \
-t 32 \
-i ont.ctg.lay.gz \
-fo Atont.ctg.fasta

quast \
	--large \
	-o Atont_wtdbg2_quast \
	-t 8 \
	Atont.ctg.fasta

makeblastdb \
       -in Atont.ctg.fasta \
       -out Atont.ctg_db \
       -dbtype nucl \
       -parse_seqids

blastn \
        -db Atont.ctg_db \
        -query AtChloRef.fasta \
        -evalue 0.1 \
        -outfmt "6 sacc evalue bitscore" \
        -out blast_out.tsv

minimap2 -ax map-ont Atont.ctg.fasta ERR2173373.fastq.gz > ontasm.sam

samtools sort -o ontasm.bam -@ 8 ontasm.sam && \
samtools index ontasm.bam && \
samtools coverage -o coverage.tsv ontasm.bam

Perl PLCL1.pl

#!/usr/bin/env perl -w
use strict;
use warnings;

my $coverage="coverage.tsv";
my $blast="blast_out.tsv";
my $outfile="blast_cov.tsv";
my $outfile2="hitcontig_cov.tsv";
my $outfile3="hitcontig_cov.u.csv";
#step1 coverage data save as hash, key is rname
my %cov_data = ();
open COV,"<coverage.tsv" or die "!";
while(<COV>){
chomp;
my($rname,$start,$end,$numreads,$covbases,$coverage,$meandepth,$meanbaseq,$meanmapq) = split (/\t/, $_, 9);
$cov_data{$rname} = [$start,$end,$numreads,$covbases,$coverage,$meandepth,$meanbaseq,$meanmapq];
}
close(COV);

#step2 merge blast data and meandepth
open OUT, '>', $outfile or die "!";
open OUT2, '>', $outfile2 or die "!";
open BLAST, '<',$blast or die "!";
while (my $line = <BLAST>){
chomp($line);
#$line =~ /^#/ and next;
my ($rname, $evalue, $bitscore, $length) = split(/\t/, $line, 4);
#search coverage ids
my $ref_data = $cov_data{$rname};
#add meandepth 1st place of cov_data
my $depth = $ref_data->[5];
#output
print OUT join ("\t",map $_ // "", $rname,$evalue,$bitscore,$length, $depth), "\n";
print OUT2 join (",", map $_ // "", $rname,$depth), "\n";
}
close(OUT);
close(OUT2);
close(BLAST);

#step3 create table with only unique value
open IN, '<'. $outfile2 or die "!";
open OUT3, '>', $outfile3 or die "!";
my @array = <IN>;
my %count;
@array=grep{!$count{$_}++}@array;
print OUT3 @array, "\n";
close(IN);
close(OUT3);

Perl PLCL2.auto.pl 2>&1 > PLCL2.auto.log

#!/usr/bin/env perl -w
use strict;
use warnings;
use Algorithm::KMeans;
#perlbrew switch perl-5.34.0
#perlbrew switch-off
#export LD_LIBRARY_PATH="/ gsl-2.6/lib"
#export LD_LIBRARY_PATH="/gsl-2.6/"
#/usr/sbin/setenforce 0

my $datafile="hitcontig_cov.u.csv";
my $mask="N1";
my $clusterer = Algorithm::KMeans -> new(datafile => $datafile,
                                        mask => $mask,
                                        K => 0,
                                        cluster_seeding => 'smart',
                                        terminal_output => 1,
                                        write_clusters_to_files => 1,
                                        );

$clusterer->read_data_from_file();

my ($clusters_hash,$cluster_centers_hash)=$clusterer->kmeans();
my $K_best=$clusterer->get_K_best();
my $QoCval=$clusterer->show_QoC_values();

my $outfile4="kmeans.out";
open OUT4, '>', $outfile4 or die "!";
print OUT4 "\ncluster\tcluster_center_value\tcontig\n";
foreach my $cluster_id (sort keys %{$clusters_hash}) {
    print OUT4 "\n$cluster_id\t@{$cluster_centers_hash->{$cluster_id}}\t@{$clusters_hash->{$cluster_id}}\n";
}
close (OUT4);

Perl PLCL2.manual.pl 2>&1 > PLCL2.manual.log

#!/usr/bin/env perl -w
use strict;
use warnings;
use Algorithm::KMeans;

my $datafile="hitcontig_cov.u.csv";
my $mask="N1";
my $clusterer = Algorithm::KMeans -> new(datafile => $datafile,
                                        mask => $mask,
                                        K => [k],
                                        cluster_seeding => 'smart',
                                        terminal_output => 1,
                                        write_clusters_to_files => 1,
                                        );

$clusterer->read_data_from_file();
my ($clusters_hash,$cluster_centers_hash)=$clusterer->kmeans();
my $K_best=$clusterer->get_K_best();
my $QoCval=$clusterer->show_QoC_values();

my $outfile4="kmeans.k[k].out";
open OUT4, '>', $outfile4 or die "!";
print OUT4 "\ncluster\tcluster_center_value\tcontig\n";
foreach my $cluster_id (sort keys %{$clusters_hash}) {
    print OUT4 "\n$cluster_id\t@{$cluster_centers_hash->{$cluster_id}}\t@{$clusters_hash->{$cluster_id}}\n";
}
close (OUT4);

#load libraries
library(ggplot2)
library(dplyr)
library(tibble)
library(ggrepel)

#load chloroplast data
chlo <- read.csv("chlo_u.km.csv")
colnames(chlo)[5] <- "kmeans_k0"
colnames(chlo)[6] <- "kmeans_k5"
chlo$coverage <- as.numeric(chlo$coverage)
colnames(chlo)

#display k=0 plot without labels
ggplot(chlo) +
  aes(x=bitscore,y=coverage)+
  geom_point(mapping = aes(color=kmeans_k0),size=4) +
  geom_point(colour="gray90",size=0.01)+
  scale_x_log10(limits=c(1,1000000),labels=scales::comma)+
  scale_y_log10(limits=c(0.1,1e4))+
  theme(text=element_text(size=20))

#display k=0 plot with labels > 50 coverage
ggplot(chlo) +
  aes(x=bitscore,y=coverage)+
  geom_point(mapping = aes(color=kmeans_k0),size=4) +
  geom_point(colour="gray90",size=0.01)+
  scale_x_log10(limits=c(1,1000000),labels=scales::comma)+
  scale_y_log10(limits=c(0.1,1e4))+
  theme(text=element_text(size=20))+
  geom_label_repel(data=chlo %>% filter(coverage>50),aes(label=tig),nudge_x=0.1,nudge_y=0.2)

#save k = 0 plot with labels > 50 coverage
c0 <- ggplot(chlo) +
  aes(x=bitscore,y=coverage)+
  geom_point(mapping = aes(color=kmeans_k0),size=4) +
  geom_point(colour="gray90",size=0.01)+
  scale_x_log10(limits=c(1,1000000),labels=scales::comma)+
  scale_y_log10(limits=c(0.1,1e4))+
  theme(text=element_text(size=20))+
  geom_label_repel(data=chlo %>% filter(coverage>50),aes(label=tig),nudge_x=0.1,nudge_y=0.2)
ggsave("At_chlo_k0.svg",c0)

#save k = 5 plot with labels > 50 coverage
c5 <- ggplot(chlo) +
  aes(x=bitscore,y=coverage)+
  geom_point(mapping = aes(color=kmeans_k5),size=4) +
  geom_point(colour="gray90",size=0.01)+
  scale_x_log10(limits=c(1,1000000),labels=scales::comma)+
  scale_y_log10(limits=c(0.1,1e4))+
  theme(text=element_text(size=20))+
  geom_label_repel(data=chlo %>% filter(coverage>50),aes(label=tig),nudge_x=0.1,nudge_y=0.2)
ggsave("At_chlo_k5.svg",c5)

mkdir $PWD/tempbam
for i in $(cat $PWD/list.txt);
do
        samtools view \
                -F 16 \
                -b \
                -o $PWD/tempbam/"$i".bam \
                $PWD/ontasm.bam \
                "$i"
done &&\
samtools merge $PWD/ontasm_filt.bam \
        -f \
        $PWD/tempbam/* && \
samtools fasta \
        $PWD/ontasm_filt.bam \
        -F 4 > \
        ontasm_filt.reads.fa

canu \
-p Atont.chlo -d Atont_chloro_canu \
usegrid=false \
genomeSize=0.15m \
minReadLength=12000 \
minOverlapLength=10000 \
corMinCoverage=8 \
trimReadsOverlap=100 \
trimReadsCoverage=100 \
rawErrorRate=0.500 \
correctedErrorRate=0.144 \
corOutCoverage=40 \
minInputCoverage=8 \
stopOnLowCoverage=8 \
-nanopore-raw ontasm_filt.reads.fa