TRITEX long-read assembly pipeline

The TRITEX pipeline needs to be run on rather powerful Unix servers. It has been tested on servers running recent CentOS Linux version controlled by the SLURM resource scheduler. Some commands like Hi-C read mapping can be submmitted as batch job. Pseudomolecule construction in R requires interactive sessions, either by direct SSH access or use of interactives SLURM sessions. In any case, the use of a terminal multiplexer such as tmux is recommended. Some jobs may run for several days. If Kerberos authentification is used, make sure that tickets are renewed.

Install the following software. The pipeline has been run successfully with the versions indicated below. Using the most recent versions should be fine as well.

The following R packages should be installed and accessible in .libPaths(). The versions specified are known to work with TRITEX. Newer versions will likely work as well.

If you want to use the R Shiny app Hi-C map inspector, you need to set up a web server running R Shiny. Refer to the R Shiny manual. Note that setting up a Shiny server is a non-trivial task of Unix systems administration and requires super-user rights.

Clone the Bitbucket repository of the TRITEX assembly pipeline. The path to the local copy of the repository will be referenced by $bitbucket in the code listings.

The scripts used below have parameters such as --bedtools or --samtools (indicated in the code listings) to specify the absolute paths to the required executables. For better readability, they are omitted from the commands.

Follow these steps if you already have a guide map. We are going to adjust it to the above format to be used in TRITEX.

If you plan on using a reference genome to build a guide map, go to the next section.

An example file is available here (MaizeB73_guidemap_marker_IBM.Rds).

You will need a BED file with the genomic positions of your markers. The fourth column should contain IDs like "seq_N", where N goes from 1 to the number of markers.

Most steps will be done in R from now on. If you are unsure whether a command should be pasted to the R or the shell prompt, hover with the mouse over the code listings. SH or R should appear in the top-right corner.

# read the BED file containing the marker regions
fread('genetic_map.bed') -> p

setnames(p, c("agp_chr", "agp_start", "agp_end", "seq"))
p[, .(css_contig=seq, popseq_cM=(agp_start+agp_end)/2/1e6, #set pseudo-cM as Mb position
    sorted_genome = as.character(NA),
    css_contig_length = agp_end - agp_start,
    popseq_alphachr = as.character(agp_chr), # it needs to be character column
    sorted_arm = NA) ] -> pp

pp[,popseq_chr:=as.integer(popseq_alphachr)]
pp[, sorted_alphachr := popseq_alphachr]
pp[, sorted_chr := popseq_chr]
pp[, sorted_lib := popseq_alphachr]

# Adjust any other columns based on your species.

# Save as RDS file for use in TRITEX pipeline
saveRDS(pp, "project_name_pseudopopseq.Rds")

In case you want to create a guide map from an available reference genome, first you need to obtain single-copy 100bp regions.

If you already have a genetic map, you can skip this step and go to the next section.

Most steps will be done in R from now on. If you are unsure whether a command should be pasted to the R or the shell prompt, hover with the mouse over the code listings. SH or R should appear in the top-right corner.

Now it is time to build the pseudo POPSEQ file (the guide map), which is used as input for the TRITEX pipeline.

Here is an example file of how your guide map should look like.

Independently of how you obtained your guide map, you should map the marker sequences to the contig assembly.

The output PAF file will be input into R into the next step.

Generate diagnostic plots to check for chimeras in the input assembly. If sequences from different chromosomes are joined, the chimeric scaffolds will have guide map markers from more than one chromosome.

# Create diagnostic plots for contigs longer than 1 Mb.
assembly$info[length >= 1e6, .(scaffold, length)][order(-length)] -> s
plot_chimeras(assembly=assembly, scaffolds=s, species="maize",
	refname="B73", autobreaks=F, mbscale=1, file="assembly_1Mb.pdf", cores=1) (1)

# Diagnostic plots for putative chimeras with strong drops in Hi-C coverage.
assembly$info[length >= 1e6 & mri <= -3, .(scaffold, length)][order(-length)] -> s
plot_chimeras(assembly=assembly, scaffolds=s, species="maize",
	refname="B73", autobreaks=F, mbscale=1, file="assembly_chimeras.pdf", cores=1)

Check the plots and search for drops in Hi-C coverage, misalignments and chimeras. For each case, assign the PDF page number from "assembly_1Mb.pdf" and the positions where it needs breaking.

For instance, in this example, the scaffolds that need to be broken are on pages 3, 6, and 9. We need to change the bins where the scaffold should be splitted.

assembly$info[length >= 1e6, .(scaffold, length)][order(-length)] -> ss

 i=3 # page number in PDF
 ss[i]$scaffold -> s
 # change the bin position observed in the plots
 assembly$cov[s, on='scaffold'][bin >= 10e6 & bin <= 22e6][order(r)][1, .(scaffold, bin)] -> b

 i=6 # page number in PDF
 ss[i]$scaffold -> s
 rbind(b, assembly$cov[s, on='scaffold'][bin >= 40e6 & bin <= 45e6][order(r)][1, .(scaffold, bin)])->b

 i=9 # page number in PDF
 ss[i]$scaffold -> s
 rbind(b, assembly$cov[s, on='scaffold'][bin >= 1e6 & bin <= 9e6][order(r)][1, .(scaffold, bin)])->b

# Rename column and plot again to double-check the breakpoints are correct
setnames(b, "bin", "br")
plot_chimeras(assembly=assembly, scaffolds=b, br=b, species="maize",
	refname="B73",  mbscale=1, file="assembly_chimeras_final.pdf", cores=1)

After checking if the break points are in the right places, proceed with breaking the scaffolds.

break_scaffolds(b, assembly, prefix="contig_corrected_v1_", slop=1e4, cores=1, species="maize") -> assembly_v2
saveRDS(assembly_v2, file="assembly_v2.Rds") (1)

Now we make the Hi-C contact matrix, or Hi-C map. Three files will be generated: a PDF with the intrachromosomal plots, a PDF with the interchromosomal plots, and a file to be loaded into R Shiny Map inspector (*_export.Rds).

fbed <- '$projectdir/assembly/project_name/project_name.p_ctg_MboI_fragments_30bp.bed'
read_fragdata(info=assembly_v2$info, file=fbed) -> frag_data

# Consider only contigs >= 1 Mb first
frag_data$info[!is.na(hic_chr) & length >= 1e6, .(scaffold, nfrag, length, chr=hic_chr, cM=popseq_cM)] -> hic_info

hic_map(info=hic_info, assembly=assembly_v2, frags=frag_data$bed,
	species="maize", ncores=1, min_nfrag_scaffold=30,
	max_cM_dist = 1000, binsize=2e5, min_nfrag_bin=10, gap_size=100) -> hic_map_v1

saveRDS(hic_map_v1, file="hic_map_v1.Rds")

# Make the Hi-C plots, files will be placed in your working directory.
snuc <- '$projectdir/assembly/project_name/project_name.p_ctg_MboI_fragments_30bp_split.nuc.txt'

hic_plots(rds="hic_map_v1.Rds", assembly=assembly_v2,
	cores=1, species="maize", nuc=snuc) -> hic_map_v2

Let’s do a second round of checking for chimeras. Broken scaffolds in the last step are renamed in the plots. It might not be necessary to break any more scaffolds; in this case, just check the plots and if everything is alright, proceed to decrease the contig size to 500 kb.

# If you'd like to  plot specific contigs:
data.table(scaffold=c('contig_corrected_v1_53', 'contig_corrected_v1_24',
	'contig_corrected_v1_5', 'contig_corrected_v1_76') -> ss

# If you'd like to check them all again:
assembly_v2$info[length >= 1e6, .(scaffold, length)][order(-length)] -> ss

plot_chimeras(assembly=assembly_v2, scaffolds=ss,
        species="maize", refname="B73", autobreaks=F,
        mbscale=1, file="assembly_v2_chimeras.pdf", cores=1)

# In case there are more chimeras, proceed like the previous step. Assign `i` with the PDF page and don't forget to change the bin. Otherwise, plot the Hi-C map with contigs >= 500 kb.

If another round of breaking scaffolds is needed, you can just run the previous block again. Don’t forget to change the output files and object names (assembly_vX, hic_map_vX).

If no more chimeras are found, proceed to decreasing the size to 500 kb.

frag_data$info[!is.na(hic_chr) & length >= 5e5, .(scaffold, nfrag, length, chr=hic_chr, cM=popseq_cM)] -> hic_info

hic_map(info=hic_info, assembly=assembly_v2, frags=frag_data$bed, species="maize",
	ncores=1, min_nfrag_scaffold=30, max_cM_dist = 1000, binsize=2e5,
	min_nfrag_bin=10, gap_size=100) -> hic_map_v2 (1)

saveRDS(hic_map_v2, file="hic_map_v2.Rds") (1)

snuc <- '$projectdir/assembly/project_name.p_ctg_MboI_fragments_30bp_split.nuc.txt'
hic_plots("hic_map_v2.Rds", assembly=assembly_v2, cores=1,
	species="maize", nuc=snuc) -> hic_map_v2 (1)


# Exporting to Excel the table that will be used for manual curation.
write_hic_map(rds="hic_map_v2.Rds", file="hic_map_v2.xlsx", species="maize")

This step is optional, yet recommended. It will help spot inverted or misplaced contigs. The collinearity plots will show marker positions against the assembly. It is possible to check the alignments and spot inverted/misplaced contigs. You can do this anytime you want to check the collinearity to the assembly. Remember that your assembly might be different compared to the reference you used, so always check the Hi-C contact maps also.

assembly_v2$cssaln[, .(scaffold, pos, guide_chr=paste0("chr", "", popseq_alphachr), guide_pos=popseq_cM)] -> x  (1)
hic_map_v2$agp[, .(scaffold, agp_start, agp_end, orientation, agp_chr)][x, on="scaffold"]-> x (2)

x[orientation == 1, agp_pos :=  agp_start + pos]
x[orientation == -1, agp_pos :=  agp_end - pos]
x[agp_chr == guide_chr] -> x

pdf("hic_map_v2_3_vs_assembly_v2.pdf", height=3*3, width=3*3)
par(mar=c(5,5,3,3))
lapply(chrNames(species="maize", agp=T)$agp_chr[1:10], function(i){   (3)
print(i)
x[i, plot(pch=".", main=i, agp_pos/1e6, guide_pos/1e6, xlab="Hi-C AGP", ylab="guide", las=1, bty='l'), on="agp_chr"]
hic_map_v2$agp[agp_chr == i & gap == T][, abline(v=agp_start/1e6, col="gray")] (2)
})
dev.off()

You can use Map Inspector Shiny app or the intrachromosomal plots for manual curation. Genomic regions containing putatively chimeric scaffolds can be clicked on to get their names and pinpoint the sites of culpable misjoins.

The source code required for deploying Hi-C Map Inspector on an R Shiny server is contained in $bitbucket/R/map_inspector.Rmd. You need to modify the paths to the data directory and the default map object.

You’ll want to check for inversions near scaffold borders (grey lines).

When modifications are spotted, change the Excel file generated in the last steps.

If an inverted contig is spotted, the column "new_orientation" on the table must be changed (Fig. 2A).

If there is a misplaced contig, the row containing it should be moved (don’t forget to change the bin number order) (Fig. 2B).

If there is an extra contig, you can remove it from the chromosome sheet and put it in the "chrUn" sheet. This will move this contig to the unanchored contig portion of the genome.

After making the changes, import the Excel file again and proceed with mapping.

read_hic_map(rds="hic_map_v2.Rds", file="hic_map_v2_edit.xlsx") -> nmap
diff_hic_map(rds="hic_map_v2.Rds", nmap, species="maize")
hic_map(species="maize", agp_only=T, map=nmap) -> hic_map_v3

saveRDS(hic_map_v3, file="hic_map_v3.Rds")

# Final Hi-C plots. You can check if the inverted scaffolds are correct.
snuc <- '$projectdir/assembly/project_name.p_ctg_MboI_fragments_30bp_split.nuc.txt'
hic_plots(rds="hic_map_v3.Rds", assembly=assembly_v2,
	cores=1, species="maize", nuc=snuc) -> hic_map_v3

If there are still wrong contigs, you can make new modifications to the Excel file and repeat this block until there are no more changes needed. You can make another collinearity plot to help you.