#! /bin/bash

########################################################################
### 		   					Assembly					   		 ###
########################################################################

## FIRST STEP: CLEANING
# trimmomatic

java -jar ~/bin/Trimmomatic-0.36/trimmomatic-0.36.jar PE -threads 3 -phred33 "$ID"_R1.fastq.gz "$ID"_R2.fastq.gz "$ID"_R1_cleaned.fastq.gz "$ID"_R1_cleaned_unpaired.fastq.gz "$ID"_R2_cleaned.fastq.gz "$ID"_R2_cleaned_unpaired.fastq.gz ILLUMINACLIP:adapters/P_up.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:50 2>"$ID".trimmomatic.log


## SECOND STEP: ASSEMBLY
# SOAPdenovo

echo "max_rd_len=150" > ./config.file
echo '[LIB]' >>./config.file
echo 'avg_ins=500' >>./config.file
echo 'asm_flag=3' >>./config.file
echo 'reverse_seq=0' >>./config.file
echo "q1="$ID"_R1_cleaned.fastq.gz" >>./config.file
echo "q2="$ID"_R2_cleaned.fastq.gz" >>./config.file
~/bin/SOAPdenovo-63mer all -d 1 -D 2 -s ./config.file -K 31 -p 10 -o "$ID" 2>SOAPdenovo.log
~/bin/GapCloser -b ./config.file -t 10 -l 150 -a "$ID".scafSeq -o "$ID".GapCloser 2>GapCloser.log

########################################################################
### 		   		      CDS identification				   		 ###
########################################################################

## FIRST STEP: BLAST 

tblastn -db $ID.GapCloser.renamed.fasta -query prot_P_Xuthus.faa -evalue 10e-20 -outfmt 6 > "$ID"_blast_out.txt

## SECOND STEP: SORT BLAST RESULTS 

Rscript ~/bin/sort_Blast.R $ID

## THIRD STEP: EXTRACT SEQUENCES FROM ASSEMBLY 

~/bin/copySeqList.py "$ID"_list.txt "$ID".GapCloser.fasta > "$ID".fasta
~/bin/renameFasta2.py "$ID" #create one file per protein .fas

## FOURTH STEP: CONSENSUS AND TRANSLATION

for line in $(ls *.fas)
do
if grep -q $ID $protein; then
prot=$line
~/bin/extract_esp.py ../$prot $ID ../Pxut_1.0_cds_sorted.fasta #create a tmp file with the reference for each contig found 
for line in $(ls tmp*)
do
~/bin/translatorx_vLocal.pl -i $line -o $line -p F -t F -c 1 #align with the reference
~/bin/removeGapUsingRefSeq "$line".nt_ali.fasta fasta Papilio_Xuthus 
done
cat *.sansGap.fst > "$prot".fasta
rm *.nt_ali.fasta
rm *.sansGap.fst
rm tmp*
~/bin/consensus.py "$prot".fasta #create a consensus with N if conflict found between contigs (example heterozygosity or sequencing error)
rm "$prot".fasta
transeq -sequence "$prot".fasta.cons.fst -outseq "$prot"_trslate.fst 
cd ../
fi
done
done

########################################################################
####                  Cross contamination cleaning                  ####
########################################################################

## FIRST STEP: FILES PREPARATION 

for line in $(ls *_trslate.fst)
do
prot=$line
~/bin/Orthofinder2.py $prot $ID # create a fasta file with the following format >"$ID"_"$prot"
done

for line in $(ls *OF*) 
do 
ID=${line%_O*} 
mv $line "$ID".fasta
done

for line in $(ls *R1*) 
do 
ID=${line%_R*} 
mv $line "$ID".L.fastq
done

for line in $(ls *R2*) 
do 
ID=${line%_R*} 
mv $line "$ID".R.fastq
done

for line in $(ls *.fasta) 
do 
sed -i -e 's/-//g' $line 
done

## SECOND STEP : CROSS CONTAMINATION WITH CROCO

~/bin/CroCo/src/CroCo_v0.1.sh --mode p --in Directory/ --threads 10 --output-prefix croco_output_Directory_ --output-level 3 --graph yes 2>&1 | tee Croco_log_file_Directory

## THIRD STEP : EXTRACT CLEAN DATAS 

cd croco_output
for line in $(cat *_list.csv)
do
echo "$line"_clean.fasta
cat "$line"_clean.fasta > "$line"_nucl.fasta
cat "$line"_lowcov.fasta >> "$line"_nucl.fasta
done

for line in $(ls *_nucl.fasta)
do
sp=${line%_nu*}
grep ">" $line > "$sp"_list_gene.csv
echo $sp >> list_sp.csv
done

for line in $(tail -n 1 list_sp.csv)
do
~/bin/select_croco_sequences.py "$line"_list_gene.csv "$line"_OF_nucl.fasta "$line"_OF_clean_nucl.fasta
done

# Note: good files are *_OF_clean_nucl.fasta with right ORF


########################################################################
####                  	  Orthology detection                  		####
########################################################################

## FIRST STEP: ORTHOFINDER

for line in $(ls *OF_clean_nucl.fasta)
do
line=${line%_nucl*}
transeq -sequence "$line"_nucl.fasta -outseq "$line"_prot.fasta # protein for OrthoFinder
done

python orthofinder.py -f Directory_with_prot_files -t 60

# SECOND STEP: EXPORT RESULTS (Proteins)

cat *OF_clean_prot.fasta > all_sp_prot.fst
~/bin/export_result_OrthoF2.py Orthogroups.txt

for line in $(ls OG*******.txt)
do
~/bin/export_fast_OrthoF2.py all_sp_prot.fst $line
done 


# SECOND STEP bis: EXPORT RESULTS (Nucleotides)

~/bin/export_result_OrthoF2.py Orthogroups.txt
cat *OF_clean_nucl.fasta > all_sp_nucl.fst

for line in $(ls OG*******.txt)
do
~/bin/export_fast_OrthoF2.py all_sp_nucl.fst $line
done 


########################################################################
####                        Alignment & rename                      ####
########################################################################

# Proteins

for line in $(ls *_prot.fst)
do
muscle -in $line -out "$line".algn -maxiters 1 -diags -sv -distance1 kbit20_3
~/bin/rename_OG.py "$line".algn ../ref_esp_tot.csv
done 

for line in $(ls listx*)
do
nohup ~/bin/align_rename.sh $line &
done

# Nucleotides

ls *_nucl.fst > list.csv
split -a 2 -n l/30 list.csv
for line in $(ls x*)
do 
mv $line list"$line"
done

for line in $(ls listx*)
do
nohup ~/bin/align_rename_nucl.sh $line &
done


########################################################################
####                          Cleaning                              ####
########################################################################

## FIRST STEP: HMMCleaner

cd Orthogroups
ls OG*| grep ".rename.fst" > list2.csv
split -a 2 -n l/30 list2.csv
for line in $(ls x*)
do 
mv $line list2"$line"
done

#prot
for gene in $(ls *.rename.fst)
do
~/bin/HMMcleanAA.pl --del-char @ $gene 5
gene=${gene%.*}
hmm=$(echo "$gene"_Hmm5.fasta)
sed -i -e  's/@//g' $hmm
sed -i -e  's/-//g' $hmm
sed -i -e  's/ /_/g' $hmm
algn=$(echo "$gene".algn.fasta)
muscle -in $hmm -out $algn -maxiters 1 -diags -sv -distance1 kbit20_3
done 


#nucl
for gene in $(ls *.rename.fst)
do
~/bin/HMMcleanNuc.pl --del-char @ $gene 5
gene=${gene%.*}
hmm=$(echo "$gene"_Hmm5.fasta)
sed -i -e  's/@//g' $hmm
sed -i -e  's/-//g' $hmm
sed -i -e  's/ /_/g' $hmm
algn=$(echo "$gene".algn.fasta)
muscle -in $hmm -out $algn -maxiters 1 -diags
done 

## SECOND STEP: TRIMAL

for line in $(ls *.rename.algn.fasta)
do
simpl=${line%.rename*}
simpl="$simpl".algn
~/bin/make_good_fasta.py $line "$simpl".fst
perl ~/bin/PaulBlock.pl "$simpl".fst 67
mv "$simpl".fst_PB67 "$simpl"_PB67.fasta
done
for line in $(ls *_PB67.fasta)
do
simpl=${line%.fasta*}
~/bin/trimAl/source/trimal -in $line -out "$simpl"_trimal.fasta -htmlout "$simpl"_trimal.fasta.html -seqoverlap 30 -resoverlap 0.75
done

for line in $(ls *_trimal.fasta)
do
~/bin/rename_trimal.py $line
done

# Note: Good files are *_prot_renamed.fasta or *_nucl_renamed.fasta


########################################################################
####        	Proteins selection for all datasets 	 	        ####
########################################################################


## DATASET 1 & 3

mkdir prot_95percent_sp
for line in $(ls *_prot_renamed.fasta)
do
~/bin/select_OrthoGroup_95percent.py $line ref_esp_tot.csv ref_esp_tot.csv tribus.csv prot_95percent_sp
done 

mkdir nucl_95percent_sp
for line in $(ls *_nucl_renamed.fasta)
do
~/bin/select_OrthoGroup_95percent.py $line ref_esp_tot.csv ref_esp_tot.csv tribus.csv nucl_95percent_sp
done 

## DATASET 2 & 4

mkdir prot_4mini_sp
for line in $(ls *_prot_renamed.fasta)
do
~/bin/select_OrthoGroup_4mini.py $line ref_esp_tot.csv ref_esp_tot.csv tribus.csv prot_4mini_sp
done 

mkdir nucl_4mini_sp
for line in $(ls *_nucl_renamed.fasta)
do
~/bin/select_OrthoGroup_4mini.py $line ref_esp_tot.csv ref_esp_tot.csv tribus.csv nucl_4mini_sp
done 


########################################################################
####       			   Matrix Concatenation	        			    ####
########################################################################


ls *_renamed.fasta > list_files.fasta
~/bin/seqCat.pl -dlist_files.fasta -if
~/bin/seqConverter.pl -dseqCat_sequences.nex -of
mv seqCat_sequences.fasta Final_matrix_Dataset1.fasta

perl ~/bin/PaulBlock.pl Final_matrix_Dataset1.fasta 99
mv Final_matrix_Dataset1.fasta_PB99 Final_matrix_Dataset1.fasta

~/bin/fasta2phylip.py Final_matrix_Dataset1.fasta > Final_matrix_Dataset1.phy


########################################################################
####            	 Supermatrix topology inference		            ####
########################################################################

## MAXIMUM LIKELIHOOD: IQ-TREE

# Dataset1_LG_BS
~/bin/iqtree -nt 64 -m LG+G4+F -b 100 -s Final_matrix_Dataset1.phy

# Dataset1_LG_UFBS
~/bin/iqtree -nt 64 -m LG+G4+F -bb 1000 -s Final_matrix_Dataset1.phy

# Dataset1_PMSF_UFBS
~/bin/iqtree -nt 18 -s Final_matrix_Dataset1.phy -m LG+C20+F+G -ft guide_tree.tre -mem 50G
~/bin/iqtree -nt 18 -s Final_matrix_Dataset1.phy -m LG+C20+F+G -fs Final_matrix_Dataset1.sitefreq -mem 50G -bb 1000

# Dataset2_LG_UFBS
~/bin/iqtree -nt 64 -m LG+G4+F -bb 1000 -s Final_matrix_Dataset2.fasta

# Dataset2_PMSF_UFBS
~/bin/iqtree -nt 64 -s Final_matrix_Dataset2.fasta -m LG+C20+F+G -ft guide_tree.tre -mem 500G -redo
~/bin/iqtree -nt 64 -s Final_matrix_Dataset2.fasta -m LG+C20+F+G -fs Final_matrix_Dataset2.fasta.sitefreq -mem 500G -bb 1000

# Dataset3_MFP_UFBS 
~/bin/iqtree -s Final_matrix_Dataset3.phy -spp RAxML_partition_genes.txt -m MFP -nt 24 -bb 1000

# Dataset4_MFP_UFBS 
~/bin/iqtree -s Final_matrix_Dataset4.phy -spp RAxML_partition_genes.txt -m MFP -nt 25 -bb 1000


## BAYESIAN: PHYLOBAYES

# Phylobayes_760_CAT_F81
mpirun -np 32 ~/bin/pbmpi-master/data/pb_mpi -d Final_matrix_Dataset1.phy -cat -dc -poisson -S -incinit 1 -x 1 10000 pbmpi1 
mpirun -np 32 ~/bin/pbmpi-master/data/pb_mpi -d Final_matrix_Dataset1.phy -cat -dc -poisson -S -incinit 1 -x 1 10000 pbmpi2 


########################################################################
####       	       	  Supertree topology inference       	    	####
########################################################################

## FIRST STEP: IQTREE

for line in $(ls *_prot_renamed.fasta)
do
~/bin/iqtree -nt 1 -m LG+G4+F -b 100 -s $line
done

## SECOND  STEP: ASTRAL

for tree in *fasta.contree; 
do 
cat $tree >> all.tre; 
echo  >> all.trees; 
done

~/bin/collapseBootstrap trees.file=all.tre bootstrap.treshold=70 # collapse branches with a bootstrap threshold < 70
java -jar astral.5.6.3.jar -t 1 -i all_collapsedAt70.txt

##  THIRD STEP: SUPERTRIPLETS

# bppReRoot is part of the bppSuite http://biopp.univ-montp2.fr/wiki/index.php/BppSuite
# outgroups.txt is a file containing two lines: "Choristoneura_fumiferana" and "Bombyx_mori"

~/bin/bppReRoot input.list.file=all_collapsedAt70.txt outgroups.file=outgroups.txt output.trees.file=all_collapsedAt70_rooted.txt  tryAgain.option=false; 
java -jar -Xmx500m  SuperTriplets_v1.1.jar all_collapsedAt70_rooted.txt superTreeST.txt
java -jar -Xmx500m  SuperTriplets_v1.1.jar all_collapsedAt70_rooted.txt superMatrixTree_ST_scores.txt superMatrixTree.tre


########################################################################
####                        	Dating	                            ####
########################################################################

## PHYLOBAYES

# Dataset 1: 59 species (ingroup + outgroup)
pb -d Final_matrix_Dataset1.phy -T tree_6621_dating_59spp.tre -cal datation_760.calib -cat -gtr -dgam 4 -bd -ugam -sb -x 1 10000 Dated_tree_59spp_CATGTR_bd_ugam_soft1
pb -d Final_matrix_Dataset1.phy -T tree_6621_dating_59spp.tre -cal datation_760.calib -cat -gtr -dgam 4 -bd -ugam -sb -x 1 10000 Dated_tree_59spp_CATGTR_bd_ugam_soft2
pb -d Final_matrix_Dataset1.phy -T tree_6621_dating_59spp.tre -cal datation_760.calib -cat -gtr -dgam 4 -bd -ugam -sb -x 1 10000 Dated_tree_59spp_CATGTR_bd_ugam_soft3

# Dataset 2: 58 species (ingroup + outgroup - Parnassius imperator removed)
pb -d Final_matrix_Dataset1.phy -T tree_6621_dating_58spp.tre -cal datation_760.calib -cat -gtr -dgam 4 -bd -ugam -sb -x 1 10000 Dated_tree_58spp_CATGTR_bd_ugam_soft1
pb -d Final_matrix_Dataset1.phy -T tree_6621_dating_58spp.tre -cal datation_760.calib -cat -gtr -dgam 4 -bd -ugam -sb -x 1 10000 Dated_tree_58spp_CATGTR_bd_ugam_soft2
pb -d Final_matrix_Dataset1.phy -T tree_6621_dating_58spp.tre -cal datation_760.calib -cat -gtr -dgam 4 -bd -ugam -sb -x 1 10000 Dated_tree_58spp_CATGTR_bd_ugam_soft3

# Dataset 3: 45 species (ingroup only)
pb -d Final_matrix_Dataset1.phy -T tree_6621_dating_45spp.tre -cal datation_760.calib -cat -gtr -dgam 4 -bd -ugam -sb -x 1 10000 Dated_tree_45spp_CATGTR_bd_ugam_soft1
pb -d Final_matrix_Dataset1.phy -T tree_6621_dating_45spp.tre -cal datation_760.calib -cat -gtr -dgam 4 -bd -ugam -sb -x 1 10000 Dated_tree_45spp_CATGTR_bd_ugam_soft2
pb -d Final_matrix_Dataset1.phy -T tree_6621_dating_45spp.tre -cal datation_760.calib -cat -gtr -dgam 4 -bd -ugam -sb -x 1 10000 Dated_tree_45spp_CATGTR_bd_ugam_soft3

# Dataset 4: 44 species (ingroup only - Parnassius imperator removed)
pb -d Final_matrix_Dataset1.phy -T tree_6621_dating_44spp.tre -cal datation_760.calib -cat -gtr -dgam 4 -bd -ugam -sb -x 1 10000 Dated_tree_44spp_CATGTR_bd_ugam_soft1
pb -d Final_matrix_Dataset1.phy -T tree_6621_dating_44spp.tre -cal datation_760.calib -cat -gtr -dgam 4 -bd -ugam -sb -x 1 10000 Dated_tree_44spp_CATGTR_bd_ugam_soft2
pb -d Final_matrix_Dataset1.phy -T tree_6621_dating_44spp.tre -cal datation_760.calib -cat -gtr -dgam 4 -bd -ugam -sb -x 1 10000 Dated_tree_44spp_CATGTR_bd_ugam_soft3
