#Importing data:Multiplexed single-end FASTQ with barcodes in sequence

qiime tools import \
  --type MultiplexedSingleEndBarcodeInSequence \
  --input-path PGM20161121.fastq.gz \
  --output-path multiplexed-seqs-3.qza


#Demultiplex the reads

qiime cutadapt demux-single \
  --i-seqs multiplexed-seqs-3.qza \
  --m-barcodes-file human_oral_3_m_map.txt \
  --m-barcodes-column BarcodeSequence \
  --p-error-rate 0 \
  --o-per-sample-sequences demultiplexed-seqs-3.qza \
  --o-untrimmed-sequences untrimmed-3.qza \
  --verbose

#visualizating trimmed-seqs.qza

qiime demux summarize \
  --i-data demultiplexed-seqs-3.qza \
  --o-visualization demultiplexed-seqs-3.qzv


#Trim adapters+primers from demulitiplexed reads

qiime cutadapt trim-single \
  --i-demultiplexed-sequences demultiplexed-seqs-3.qza \
  --p-front GATGTGCCAGCMGCCGCGGTAA \
  --p-error-rate 0 \
  --p-minimum-length 200 \
  --o-trimmed-sequences trimmed-seqs-3.qza \
  --verbose


#visualizating trimmed-seqs.qza

qiime demux summarize \
  --i-data trimmed-seqs-3.qza \
  --o-visualization trimmed-seqs-3.qzv


#Sequence quality control and feature table construction

nohup qiime dada2 denoise-single \
  --p-n-threads 16 \
  --i-demultiplexed-seqs trimmed-seqs-3.qza \
  --p-trim-left 0 \
  --p-trunc-len 233 \
  --o-representative-sequences rep-seqs-3.qza \
  --o-table table-3.qza \
  --o-denoising-stats stats-3.qza

#visualizating stats-3.qza

qiime metadata tabulate \
  --m-input-file stats-3.qza \
  --o-visualization stats-3.qzv

# merge data ########################################################
#merge feature-table

qiime feature-table merge \
  --i-tables table-1.qza \
  --i-tables table-2.qza \
  --i-tables table-3.qza \
  --o-merged-table table.qza

#merge-seqs feature-table

qiime feature-table merge-seqs \
  --i-data rep-seqs-1.qza \
  --i-data rep-seqs-2.qza \
  --i-data rep-seqs-3.qza \
  --o-merged-data rep-seqs.qza


#feature-table summarize

qiime feature-table summarize \
  --i-table table.qza \
  --o-visualization table.qzv \
  --m-sample-metadata-file merge_map.txt

#Representative sequences

qiime feature-table tabulate-seqs \
--i-data rep-seqs.qza \
--o-visualization rep_seqs.qzv

#Filter features that only appear in a single sample

qiime feature-table filter-features \
--i-table table.qza \
--p-min-samples 2 \
--o-filtered-table table1.qza

qiime feature-table summarize \
  --i-table table1.qza \
  --o-visualization table1.qzv 


#Filter features by reads that appears 10 times or fewer across all samples

qiime feature-table filter-features \
--i-table table1.qza \
--p-min-frequency 10 \
--o-filtered-table table2.qza

qiime feature-table summarize \
  --i-table table2.qza \
  --o-visualization table2.qzv 

 
  
# Taxonomic analysis ########################################################
########### Taxonomic analysis---SILVA_132_QIIME_release#########
mkdir taxonomy
cd taxonomy

#assign taxonomy to features

mkdir SILVA
cd SILVA

#Get SILVA database

wget https://www.arb-silva.de/fileadmin/silva_databases/qiime/Silva_132_release.zip
unzip Silva_132_release.zip
rm Silva_132_release.zip


#Read database into QIIME2

cd SILVA_132_QIIME_release

#The sequences

qiime tools import \
--type 'FeatureData[Sequence]' \
--input-path rep_set/rep_set_16S_only/99/silva_132_99_16S.fna \
--output-path 99_otus_16S.qza

#The taxonomy strings

qiime tools import \
--type 'FeatureData[Taxonomy]' \
--input-format HeaderlessTSVTaxonomyFormat \
--input-path taxonomy/16S_only/99/majority_taxonomy_7_levels \
--output-path 99_otus_16S_taxonomy.qza

#Extract the V3/V4 region from the reference database
#--p-f-primer the primer on the left, 515F
#--p-r-primer the primer on the right, 806R
#--p-trunc-len  this value must be equal or less than the length of your trimed data

qiime feature-classifier extract-reads \
--i-sequences 99_otus_16S.qza \
--p-f-primer GTGCCAGCMGCCGCGGTAA \
--p-r-primer GGACTACHVGGGTWTCTAAT \
--p-trunc-len 233 \
--p-min-length 200 \
--p-max-length 250 \
--o-reads ref_seqs_16S_V3-V4-233.qza 

#this produces the QIIME2 artifact ref_seqs_16S_V3-V4.qza; the V3/V4 reference sequences.
#Train the classifier on this region

qiime feature-classifier fit-classifier-naive-bayes \
--i-reference-reads ref_seqs_16S_V3-V4-233.qza \
--i-reference-taxonomy 99_otus_16S_taxonomy.qza \
--o-classifier classifier_16S_V3-V4-233.qza 

#This step took a couple of hours. This produces the QIIME2 artifact classified_rep_seqs.qza.

cd taxonomy_silva132

#Classify rep seqs

qiime feature-classifier classify-sklearn \
--i-classifier taxonomy/classifier_16S_V3-V4-233.qza \
--i-reads merge_data/rep-seqs.qza \
--o-classification taxonomy/classified_rep_seqs-233.qza

#Tabulate the features, their taxonomy and the confidence of taxonomy assignment

qiime metadata tabulate \
--m-input-file classified_rep_seqs-233.qza \
--o-visualization classified_rep_seqs-233.qzv

#classified_rep_seqs barplot---table2

qiime taxa barplot \
  --i-table merge_data/table2.qza \
  --i-taxonomy taxonomy/classified_rep_seqs-233.qza \
  --m-metadata-file merge_data/merge_map.txt \
  --o-visualization taxonomy/taxabarplots_table2.qzv


############# Taxonomic analysis---gg-13-8-99-515-806-nb-classifier ############
qiime feature-classifier classify-sklearn \
  --i-classifier gg-13-8-99-515-806-nb-classifier.qza \
  --i-reads rep-seqs.qza \
  --o-classification taxonomy.qza


############# Taxonomic analysis---silva-132-99-515-806-nb-classifier ############

qiime feature-classifier classify-sklearn \
  --i-classifier silva-132-99-515-806-nb-classifier.qza \
  --i-reads rep-seqs.qza \
  --o-classification taxonomy.qza

#classified_rep_seqs barplot

qiime metadata tabulate \
  --m-input-file taxonomy.qza \
  --o-visualization taxonomy.qzv

qiime taxa barplot \
  --i-table table.qza \
  --i-taxonomy taxonomy.qza \
  --m-metadata-file merge_map.txt \
  --o-visualization taxa-bar-plots.qzv

mkdir collapse

#Collapse table at phylum level

qiime taxa collapse \
--i-table table2.qza \
--i-taxonomy classified_rep_seqs-233.qza \
--p-level 2 \
--o-collapsed-table collapse/feature_table_level2.qza

#Collapse table at class level

qiime taxa collapse \
--i-table table2.qza \
--i-taxonomy classified_rep_seqs-233.qza \
--p-level 3 \
--o-collapsed-table collapse/feature_table_level3.qza

#Collapse table at order level

qiime taxa collapse \
--i-table table2.qza \
--i-taxonomy classified_rep_seqs-233.qza \
--p-level 4 \
--o-collapsed-table collapse/feature_table_level4.qza

#Collapse table at family level

qiime taxa collapse \
--i-table table2.qza \
--i-taxonomy classified_rep_seqs-233.qza \
--p-level 5 \
--o-collapsed-table collapse/feature_table_level5.qza

#Collapse table at genus level

qiime taxa collapse \
--i-table table2.qza \
--i-taxonomy classified_rep_seqs-233.qza \
--p-level 6 \
--o-collapsed-table collapse/feature_table_level6.qza
 
 
qiime tools extract \
--input-path collapse/feature_table_level6.qza \
--output-path collapse/feature_table_level6


biom convert \
-i collapse/feature_table_level6/4e6797a1-cb88-472e-a906-798e8ff496ab/data/feature-table.biom \
-o lcollapse/feature_table_level6.tsv \
--to-tsv

#Generate a tree for phylogenetic diversity analyses

qiime phylogeny align-to-tree-mafft-fasttree \
  --p-n-threads 16 \
  --i-sequences rep-seqs.qza \
  --o-alignment aligned-rep-seqs.qza \
  --o-masked-alignment masked-aligned-rep-seqs.qza \
  --o-tree unrooted-tree.qza \
  --o-rooted-tree rooted-tree.qza
  
#Alpha and beta diversity analysis---table2

qiime diversity core-metrics-phylogenetic \
  --i-phylogeny rooted-tree.qza \
  --i-table table2.qza \
  --p-sampling-depth 10738 \
  --m-metadata-file merge_map.txt \
  --output-dir core-metrics-results-table2

#Alpha diversity analysis---faith

qiime diversity alpha-group-significance \
  --i-alpha-diversity core-metrics-results-table2/faith_pd_vector.qza \
  --m-metadata-file merge_map.txt \
  --o-visualization core-metrics-results-table2/faith-pd-group-significance.qzv

#Alpha diversity analysis---eveness

qiime diversity alpha-group-significance \
  --i-alpha-diversity core-metrics-results-table2/evenness_vector.qza \
  --m-metadata-file merge_map.txt \
  --o-visualization core-metrics-results-table2/evenness-group-significance.qzv


#Alpha diversity analysis---shannon

qiime diversity alpha-group-significance \
  --i-alpha-diversity core-metrics-results-table2/shannon_vector.qza \
  --m-metadata-file merge_map.txt \
  --o-visualization core-metrics-results-table2/shannon_vector.qzv
  
#Alpha diversity analysis---observed_otus

qiime diversity alpha-group-significance \
  --i-alpha-diversity core-metrics-results-table2/observed_otus_vector.qza \
  --m-metadata-file merge_map.txt \
  --o-visualization core-metrics-results-table2/observed_otus_vector.qzv




#Beta diversity analysis---unweighted_unifrac_distance_matrix

qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results-table2/unweighted_unifrac_distance_matrix.qza \
  --m-metadata-file merge_map.txt \
  --m-metadata-column Treatment \
  --o-visualization core-metrics-results-table2/unweighted-unifrac-Treatment-significance.qzv \
  --p-pairwise

#Beta diversity analysis---jaccard_distance_matrix

qiime diversity beta-group-significance \
  --i-distance-matrix jaccard_distance_matrix.qza \
  --m-metadata-file merge_map.txt \
  --m-metadata-column Treatment \
  --o-visualization jaccard-Treatment-significance.qzv \
  --p-pairwise
  
#Beta diversity analysis---bray_curtis_distance_matrix

qiime diversity beta-group-significance \
  --i-distance-matrix bray_curtis_distance_matrix.qza \
  --m-metadata-file merge_map.txt \
  --m-metadata-column Treatment \
  --o-visualization bray_curtis-Treatment_significance.qzv \
  --p-pairwise

#Beta diversity analysis---weighted_unifrac_distance_matrix

qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics-results-table2/weighted_unifrac_distance_matrix.qza \
  --m-metadata-file merge_map.txt \
  --m-metadata-column Treatment \
  --o-visualization core-metrics-results-table2/weighted-unifrac-Treatment-significance.qzv \
  --p-pairwise



#emperor tool: explore principal coordinates (PCoA) plots

qiime emperor plot \
  --i-pcoa core-metrics-results-table2/unweighted_unifrac_pcoa_results.qza \
  --m-metadata-file merge_map.txt \
  --o-visualization core-metrics-results-table2/unweighted-unifrac_pcoa_results.qzv


qiime emperor plot \
  --i-pcoa core-metrics-results-table2/bray_curtis_pcoa_results.qza \
  --m-metadata-file merge_map.txt \
  --o-visualization core-metrics-results-table2/bray-curtis_pcoa_results.qzv




#Alpha rarefaction plotting---table2

qiime diversity alpha-rarefaction \
  --i-table table2.qza \
  --i-phylogeny rooted-tree.qza \
  --p-max-depth 51000\
  --m-metadata-file merge_map.txt \
  --o-visualization alpha-rarefaction-table2-51000.qzv

# ANCOM ##########
mkdir ANCOM

#ANCOM for original otu_table---table2

qiime composition add-pseudocount \
  --i-table table2.qza \
  --o-composition-table ANCOM/comp-table2.qza

#Run ANCOM

qiime composition ancom \
  --i-table ANCOM/comp-table2.qza \
  --m-metadata-file merge_map.txt \
  --m-metadata-column Treatment \
  --o-visualization ANCOM/ancom-Treatment-table2.qzv


#Differential abundance testing with ANCOM --- Treatment="ND"~Sex

qiime feature-table filter-samples \
  --i-table table2.qza \
  --m-metadata-file merge_map.txt \
  --p-where "[Treatment]='ND'" \
  --o-filtered-table ANCOM/ND-table.qza

qiime composition add-pseudocount \
  --i-table ANCOM/ND-table.qza \
  --o-composition-table ANCOM/comp-ND-table.qza

nohup qiime composition ancom \
  --i-table ANCOM/comp-ND-table.qza \
  --m-metadata-file merge_map.txt  \
  --m-metadata-column Description \
  --o-visualization ANCOM/ancom-ND-Description.qzv


#Differential abundance testing with ANCOM --- Treatment="D"~Sex

qiime feature-table filter-samples \
  --i-table table2.qza \
  --m-metadata-file merge_map.txt \
  --p-where "[Treatment]='D'" \
  --o-filtered-table ANCOM/D-table.qza

qiime composition add-pseudocount \
  --i-table ANCOM/D-table.qza \
  --o-composition-table ANCOM/comp-D-table.qza

qiime composition ancom \
  --i-table ANCOM/comp-D-table.qza \
  --m-metadata-file merge_map.txt  \
  --m-metadata-column Description \
  --o-visualization ANCOM/ancom-D-Description.qzv

#ANCOM for collapse table 

qiime composition add-pseudocount \
--i-table collapse/feature_table_level6.qza \
--o-composition-table ANNCOM/ANCOM_ready_table_level6.qza

qiime composition ancom \
  --i-table ANCOM/ANCOM_ready_table_level6 \
  --m-metadata-file merge_map.txt  \
  --m-metadata-column Treatment \
  --o-visualization ANCOM/ancom-Treatment-level6.qzv



# gneiss ##########
mkdir gneiss

#Also check remaining read counts

qiime feature-table summarize \
--i-table table2.qza \
--o-visualization gneiss/feature_table.qzv

#Correlation-clustering

qiime gneiss correlation-clustering \
  --i-table table2.qza \
  --o-clustering gneiss/hierarchy.qza

#computes the log ratios between groups at each node in the tree.

qiime gneiss ilr-hierarchical \
  --i-table table2.qza \
  --i-tree gneiss/hierarchy.qza \
  --o-balances gneiss/balances.qza

#run linear regression on the balances

qiime gneiss ols-regression \
  --p-formula "Treatment" \
  --i-table gneiss/balances.qza \
  --i-tree gneiss/hierarchy.qza \
  --m-metadata-file merge_map.txt \
  --o-visualization gneiss/regression_summary_Treatment.qzv

#run linear regression on the balances---different factors
qiime gneiss ols-regression \
  --p-formula "Diabetes+Age+HLA_Risk_Group" \
  --i-table gneiss2/balances.qza \
  --i-tree gneiss2/hierarchy.qza \
  --m-metadata-file new_metadata.txt \
  --o-visualization gneiss2/regression_summary_threefactors.qzv

qiime gneiss ols-regression \
  --p-formula "Description" \
  --i-table gneiss/balances.qza \
  --i-tree gneiss/hierarchy.qza \
  --m-metadata-file merge_map.txt \
  --o-visualization gneiss/regression_summary_Description.qzv

#visualize these balances on a heatmap

qiime gneiss dendrogram-heatmap \
  --i-table table2.qza \
  --i-tree gneiss/hierarchy.qza \
  --m-metadata-file merge_map.txt \
  --m-metadata-column Treatment \
  --p-color-map seismic \
  --o-visualization gneiss/heatmap.qzv

#plot a boxplot and identify taxa that could be explaining the differences between the D and ND groups.

qiime gneiss balance-taxonomy \
  --i-table table2.qza \
  --i-tree gneiss/hierarchy.qza \
  --i-taxonomy classified_rep_seqs-233.qza \
  --p-taxa-level 2 \
  --p-balance-name 'y0' \
  --m-metadata-file merge_map.txt \
  --m-metadata-column Treatment \
  --o-visualization gneiss/y0_taxa_summary_lv2.qzv

nohup qiime gneiss balance-taxonomy \
  --i-table table2.qza \
  --i-tree gneiss/hierarchy.qza \
  --i-taxonomy classified_rep_seqs-233.qza \
  --p-taxa-level 6 \
  --p-balance-name 'y0' \
  --m-metadata-file merge_map.txt \
  --m-metadata-column Treatment \
   --o-visualization gneiss/y0_taxa_summary_lv6.qzv

qiime tools extract \
--input-path rooted_tree.qza \
--output-path rooted_tree

qiime tools extract \
--input-path classified_rep_seqs-233.qza \
--output-path classified_rep_seqs-233

qiime tools extract \
--input-path table.qza \
--output-path table

biom convert \
-i table2/9c40c08e-45c1-405c-bc41-27c3e05025c5/data/feature-table.biom \
-o table2/out_table.tsv \
--to-tsv

biom convert -h

# small tips  ########################################################
#Count the number of sequences in a FASTQ file

grep -n filename.fastq

#Count the total number of rows in a file that contain primers

grep -c ‘primer/barcode/adapter/’ filename.fastq

#count the total number of rows in a file

wc -l filename.fastq

#compress file

gzip filename
gzip -c filename

#unzip file

gzip -d filename
gunzip filename

#check .fastq file

less -S -N filename.fastq
zless -S -N filename.fastq.gz

cat -n filename.fastq
zcat -n filename.fastq.gz
zcat filename.fastq.gz | head -n 10

head -n filename.fastq
tail -n filename.fastq




biom add-metadata \
-i exported/feature-table.biom \
-o table-with-taxonomy.biom \
--observation-metadata-fp exported/taxonomy.tsv \
--observation-header taxonomy \
--sc-separated taxonomy

biom convert \
-i table-with-taxonomy.biom \
-o table-with-taxonomy.tsv \
--to-tsv \
--header-key taxonomy


biom add-metadata \
-i feature-table.biom \
-o table- with -taxonomy.biom \
–-observation-metadata-fp biom-taxonomy.tsv \
--sc-separated taxonomy