4.2.1 Run BLAST on server

Using a qsub run the following code (takes about 2 days - breaks in the middle because takes too much time)We use the OTU file generated at step 06.

DIR_PROJECT="/projet/sbr/ccebarcodep1408/nz_2017/blast/"

cd $DIR_PROJECT

FILE=nz_2017_otu_mothur

FASTA=$DIR_PROJECT$FILE".fas"
BLAST_XLM=$DIR_PROJECT$FILE".xml"
BLAST_TSV=$DIR_PROJECT$FILE".tsv"
OUT_FMT="6 qseqid sseqid sacc stitle sscinames staxids sskingdoms sblastnames pident slen length mismatch gapopen qstart qend sstart send evalue bitscore"

# blastall -p blastn -b 100 -v 100 -m 7 -e 1.00e-100 -i $FASTA -o $BLAST_XLM -d /db/blast/all/nt
blastn -max_target_seqs 100 -evalue 1.00e-10 -query $FASTA -out $BLAST_TSV -db /db/blast/all/nt -outfmt "$OUT_FMT"

4.2.2 Process BLAST output

We call the blast_18S function to create two files

• a modified BLAST output file which includes additional coolumns
  • kingdom -> species : PR2 taxonomy for those accession numbers that are present in PR2
  • hit_rank : the rank of the hit based on decreasing % identity and decreasing bit scores
  • uncultured : TRUE if the hit corresponds to an uncultured item
  • hit_lineage : GenBank taxonomy of the hit
• a summary file with one line per query that contains three sets of columns
  • The top hit (column with prefix hit_top_)
  • The top hit for which a PR2 sequence is available (columns starting with hit_pr2_).
  • The top hit corresponding to a culture or an isolate (columns starting with hit_cult_).
  • A “consensus”" for taxonomy at the genus level based on sequences which have > 96% similairty and for which PR2 annotation exist.
  • A remark concerning any contradiction at the division level between hits that have > 90% similarity. This may help detect chimera if the balance is 50:50. If this only for a small number of hits, this no problem (bad annotation ?)

blast_file_name <- str_c(path_blast, "nz_2017_otu_mothur.blast.tsv")
blast_18S_reformat(blast_file_name)

Read the modified BLAST file (100 lines per otu) and the BLAST summary file (1 line per otu).

blast_file <- str_c(path_blast, "nz_2017_otu_mothur.blast_pr2.tsv")
blast <- read_tsv(blast_file, guess_max = 10000)
colnames(blast)

 [1] "query_id"           "hit_id"             "hit_acc"           
 [4] "hit_title"          "hit_sci_names"      "hit_tax_ids"       
 [7] "hit_super_kingdoms" "hit_blast_names"    "pct_identity"      
[10] "hit_length"         "alignment_length"   "mismatches"        
[13] "gap_opens"          "query_start"        "query_end"         
[16] "hit_start"          "hit_end"            "evalue"            
[19] "bit_score"          "hit_rank"           "species"           
[22] "kingdom"            "supergroup"         "division"          
[25] "class"              "order"              "family"            
[28] "genus"              "uncultured"         "tax_organism"      
[31] "tax_lineage"       

blast_summary_file <- str_c(path_blast, "nz_2017_otu_mothur.blast_summary.tsv")
blast_summary <- read_tsv(blast_summary_file, guess_max = 10000)
colnames(blast_summary)

 [1] "query_id"                    "hit_number"                 
 [3] "hit_top_hit_id"              "hit_top_hit_acc"            
 [5] "hit_top_hit_title"           "hit_top_hit_sci_names"      
 [7] "hit_top_hit_tax_ids"         "hit_top_hit_super_kingdoms" 
 [9] "hit_top_hit_blast_names"     "hit_top_pct_identity"       
[11] "hit_top_hit_length"          "hit_top_alignment_length"   
[13] "hit_top_mismatches"          "hit_top_gap_opens"          
[15] "hit_top_query_start"         "hit_top_query_end"          
[17] "hit_top_hit_start"           "hit_top_hit_end"            
[19] "hit_top_evalue"              "hit_top_bit_score"          
[21] "hit_top_hit_rank"            "hit_top_species"            
[23] "hit_top_kingdom"             "hit_top_supergroup"         
[25] "hit_top_division"            "hit_top_class"              
[27] "hit_top_order"               "hit_top_family"             
[29] "hit_top_genus"               "hit_top_uncultured"         
[31] "hit_top_tax_organism"        "hit_top_tax_lineage"        
[33] "hit_pr2_hit_id"              "hit_pr2_hit_acc"            
[35] "hit_pr2_hit_title"           "hit_pr2_hit_sci_names"      
[37] "hit_pr2_hit_tax_ids"         "hit_pr2_hit_super_kingdoms" 
[39] "hit_pr2_hit_blast_names"     "hit_pr2_pct_identity"       
[41] "hit_pr2_hit_length"          "hit_pr2_alignment_length"   
[43] "hit_pr2_mismatches"          "hit_pr2_gap_opens"          
[45] "hit_pr2_query_start"         "hit_pr2_query_end"          
[47] "hit_pr2_hit_start"           "hit_pr2_hit_end"            
[49] "hit_pr2_evalue"              "hit_pr2_bit_score"          
[51] "hit_pr2_hit_rank"            "hit_pr2_species"            
[53] "hit_pr2_kingdom"             "hit_pr2_supergroup"         
[55] "hit_pr2_division"            "hit_pr2_class"              
[57] "hit_pr2_order"               "hit_pr2_family"             
[59] "hit_pr2_genus"               "hit_pr2_uncultured"         
[61] "hit_pr2_tax_organism"        "hit_pr2_tax_lineage"        
[63] "hit_cult_hit_id"             "hit_cult_hit_acc"           
[65] "hit_cult_hit_title"          "hit_cult_hit_sci_names"     
[67] "hit_cult_hit_tax_ids"        "hit_cult_hit_super_kingdoms"
[69] "hit_cult_hit_blast_names"    "hit_cult_pct_identity"      
[71] "hit_cult_hit_length"         "hit_cult_alignment_length"  
[73] "hit_cult_mismatches"         "hit_cult_gap_opens"         
[75] "hit_cult_query_start"       
 [ reached getOption("max.print") -- omitted 27 entries ]

• Split the query_id column into the otuxxxx field and the taxonomy.
  • Note | is a metacharacter, so it must be escaped.
• Save a file with the current taxonomy. This file will be used to correct the taxonomy

mothur_otus <- str_split_fixed(blast_summary$query_id, pattern = "\\|", 8)
colnames(mothur_otus) <- c("OTUNumber", str_c("mothur_", pr2_taxo_levels[2:8]))
blast_summary <- cbind(mothur_otus, blast_summary)
write_tsv(as.data.frame(mothur_otus), str_c(path_mothur, "nz_2017_otu_mothur_taxo.tsv"), 
    na = "")

4.3.1 Percent identity for best hit

• Mean identity of best hit (%) = 98.11

plot_pct <- gg_hist(blast_summary, "hit_top_pct_identity")
plot_pct + xlim(100, 80) + ylim(0, 0.75) + ggtitle("Best hit")

summary(blast_summary$hit_top_pct_identity, digits = 4)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  77.66   97.60   99.40   98.11  100.00  100.00 

• Mean identity of best hit to isolates (%) = 93.80

plot_pct <- gg_hist(blast_summary, "hit_cult_pct_identity")
plot_pct + xlim(100, 80) + ylim(0, 0.25) + ggtitle("Best hit to isolate")

summary(blast_summary$hit_cult_pct_identity, digits = 4)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
  74.24   90.42   94.72   93.80   98.15  100.00    1422 

4.3.2 Average percent identity vs OTU rank

• The percent identity becomes lower for the less abundant OTUS

ggplot(aes(x = as.numeric(row.names(blast_summary)), y = hit_top_pct_identity), 
    data = blast_summary) + geom_point() + geom_smooth() + ylim(100, 80) + labs(x = "otu number", 
    y = "Top hit identity (%)")

4.3.3 Average percent identity for different divisions

The idea is to compare the % identity for the different divisions, with the idea that for some divisions we have very few reference sequences and therefore probably new sequences. This done for the first 100 OTUs and for all OTUs

ggplot(aes(x = mothur_division, y = hit_top_pct_identity), data = blast_summary) + 
    geom_boxplot() + ylim(100, 80) + coord_flip() + labs(x = "Division", y = "Top hit identity (%)") + 
    ggtitle("For all Otus")

ggplot(aes(x = mothur_division, y = hit_top_pct_identity), data = blast_summary[1:100, 
    ]) + geom_boxplot() + ylim(100, 96) + coord_flip() + labs(x = "Division", 
    y = "Top hit identity (%)") + ggtitle("For first 100 Otus")

4.5.1 Merge the dada2 assignement

Because of this problem we use the dada2 assignement rather than the mothur assignement to do the comparison with the BLAST.

• First we merge the dada2 assignement (see Step 06)

phyloseq_dada2 <- readRDS(phyloseq, file = "phyloseq_nz_2017_dada2.rds")

taxo_dada2 <- data.frame(tax_table(phyloseq_dada2)) %>% rename_all(funs(str_c("dada2_", 
    .))) %>% rownames_to_column(var = "OTUNumber")
blast_summary <- left_join(blast_summary, taxo_dada2)

4.5.2 Compare the dada2 taxonomy and the BLAST annotation at the genus and division level

• We create a file (NZ_2017_dada2_blast_disagree.tsv) that lists the otus for which there is a disagrement at the genus level between the dada2 taxonomy and the BLAST taxonomy

blast_disagree_genus <- blast_summary %>% filter(dada2_genus != hit_pr2_genus)
write_tsv(blast_disagree_genus, str_c(path_blast, "NZ_2017_dada2_blast_disagree.tsv"), 
    na = "")

blast_disagree_division <- blast_summary %>% filter(dada2_division != hit_pr2_division)

• Number of otus with difference in dada2 annotation and BLAST annotation (taxonomy is based on pr2). Note that this is much better than with the mothur assingation. At the division levels, more than 10 times less disagreement. At the genus level almost twice less.

  • Division level : 155

  • Genus level : 1148

4.6.1 Examples of disagreements and change

In contrast to mothur, dada2 has a tendency to attribute an otu to a species rather than an uncultured group. Much of the time it is fine except in a few cases where PR2 annotation is not accurate and an environmental sequences is attributed wrongly to a species (often the case with dinos)

• otu0029 : dada2 labelled as Pycnococcacae while it hits Pycnonoccus with 100% similarity

blast_summary_display("Otu0029")

blast_display("Otu0029")

• otu0039 : dada2 labelled as Gephyrocapsa but Emiliania huxleyi has same sequence

blast_summary_display("Otu0039")

blast_display("Otu0039")

• otu0037 : For some groups V4 is not discrimant like for clade VII A2 vs A3 so the hits are from both

blast_summary_display("Otu0037")

blast_display("Otu0037")

• otu0082 : dada2 assign as Apocalathium but the hit is quite low (90%) so better to assign to Dinophyceae_XXX_sp.

blast_summary_display("Otu0082")

blast_display("Otu0082")

4.7.1 List of changes made

taxo_corrected <- read_xlsx(str_c(path_files, "NZ_2017_BLAST.xlsx"), sheet = "taxonomy corrected")
taxo_corrected_show <- taxo_corrected %>% filter(!is.na(corrected_species)) %>% 
    select(OTUNumber, dada2_species, corrected_species, corrected_class:corrected_genus, 
        corrected_remark)

knitr::kable(head(taxo_corrected_show, 31), caption = "Taxonomic corrections made") %>% 
    kable_styling(latex_options = "scale_down", font_size = 10)

4.7.2 Correct the phyloseq file

The file is saved into phyloseq_nz_2017.rds

taxo_changed <- taxo_corrected %>% filter(!is.na(corrected_species)) %>% select(OTUNumber, 
    corrected_supergroup:corrected_species) %>% rename_all(funs(str_replace(., 
    "corrected_", "")))

taxo_kept <- taxo_corrected %>% filter(is.na(corrected_species)) %>% select(OTUNumber, 
    dada2_supergroup:dada2_species) %>% rename_all(funs(str_replace(., "dada2_", 
    "")))

taxo_final <- rbind(taxo_kept, taxo_changed) %>% arrange(OTUNumber) %>% mutate(kingdom = "Eukaryota") %>% 
    column_to_rownames(var = "OTUNumber")

phyloseq_final <- phyloseq_dada2
tax_table(phyloseq_final) <- as.matrix(taxo_final)
saveRDS(phyloseq_final, file = "phyloseq_nz_2017.rds")