#36 push all new changes and files

Test_cdhit.sh

@@ -0,0 +1,238 @@
+#!/bin/bash
+
+#SBATCH --mem 18G
+#SBATCH --job-name cd-hit-test
+#SBATCH --mail-user [email protected] ## CHANGE THIS TO YOUR EMAIL
+#SBATCH --mail-type ALL
+#SBATCH -n 2
+#SBATCH -N 1
+#SBATCH -A h3abionet
+#SBATCH -o /home/groups/h3abionet/RefGraph/results/NeginV_Test_Summer2021/slurm_output/slurm-%A.out
+
+### This tests various parameters for cdhit at the annotation step UIUC pipeline 
+## Date File Created: April 23, 2022
+
+
+# Set working directory -------
+cd /home/groups/h3abionet/RefGraph/results/NeginV_Test_Summer2021/results/annotation
+
+# Load nextflow ------
+module load CD-HIT/4.8.1-IGB-gcc-8.2.0
+
+
+# Use cd-hit to cluster and remove redundancy ------
+
+# GRCH38_decoys -----
+
+# megahit -----
+
+# 0.9
+cd-hit-est \
+-i Merged_Reads/megahit/merged_sequences_GRCH38_decoys.fasta \
+-o Cluster_CDHIT/megahit/clustered_GRCH38_decoys_n5_i0.9.fasta \
+-c 0.9 \
+-n 5 \
+-T $SLURM_NPROCS
+
+# 0.92
+cd-hit-est \
+-i Merged_Reads/megahit/merged_sequences_GRCH38_decoys.fasta \
+-o Cluster_CDHIT/megahit/clustered_GRCH38_decoys_n5_i0.92.fasta \
+-c 0.92 \
+-n 5 \
+-T $SLURM_NPROCS
+
+# 0.94
+cd-hit-est \
+-i Merged_Reads/megahit/merged_sequences_GRCH38_decoys.fasta \
+-o Cluster_CDHIT/megahit/clustered_GRCH38_decoys_n5_i0.94.fasta \
+-c 0.94 \
+-n 5 \
+-T $SLURM_NPROCS
+
+# 0.96
+cd-hit-est \
+-i Merged_Reads/megahit/merged_sequences_GRCH38_decoys.fasta \
+-o Cluster_CDHIT/megahit/clustered_GRCH38_decoys_n5_i0.96.fasta \
+-c 0.96 \
+-n 5 \
+-T $SLURM_NPROCS
+
+
+# masurca -----
+
+# 0.9
+cd-hit-est \
+-i Merged_Reads/masurca/merged_sequences_GRCH38_decoys.fasta \
+-o Cluster_CDHIT/masurca/clustered_GRCH38_decoys_n5_i0.9.fasta \
+-c 0.9 \
+-n 5 \
+-T $SLURM_NPROCS
+
+# 0.92
+cd-hit-est \
+-i Merged_Reads/masurca/merged_sequences_GRCH38_decoys.fasta \
+-o Cluster_CDHIT/masurca/clustered_GRCH38_decoys_n5_i0.92.fasta \
+-c 0.92 \
+-n 5 \
+-T $SLURM_NPROCS
+
+# 0.94
+cd-hit-est \
+-i Merged_Reads/masurca/merged_sequences_GRCH38_decoys.fasta \
+-o Cluster_CDHIT/masurca/clustered_GRCH38_decoys_n5_i0.94.fasta \
+-c 0.94 \
+-n 5 \
+-T $SLURM_NPROCS
+
+# 0.96
+cd-hit-est \
+-i Merged_Reads/masurca/merged_sequences_GRCH38_decoys.fasta \
+-o Cluster_CDHIT/masurca/clustered_GRCH38_decoys_n5_i0.96.fasta \
+-c 0.96 \
+-n 5 \
+-T $SLURM_NPROCS
+
+
+# GRCH38_p0 -----
+
+# megahit -----
+
+# 0.9
+cd-hit-est \
+-i Merged_Reads/megahit/merged_sequences_GRCH38_p0.fasta \
+-o Cluster_CDHIT/megahit/clustered_GRCH38_p0_n5_i0.9.fasta \
+-c 0.9 \
+-n 5 \
+-T $SLURM_NPROCS  
+
+# 0.92
+cd-hit-est \
+-i Merged_Reads/megahit/merged_sequences_GRCH38_p0.fasta \
+-o Cluster_CDHIT/megahit/clustered_GRCH38_p0_n5_i0.92.fasta \
+-c 0.92 \
+-n 5 \
+-T $SLURM_NPROCS
+
+# 0.94
+cd-hit-est \
+-i Merged_Reads/megahit/merged_sequences_GRCH38_p0.fasta \
+-o Cluster_CDHIT/megahit/clustered_GRCH38_p0_n5_i0.94.fasta \
+-c 0.94 \
+-n 5 \
+-T $SLURM_NPROCS  
+
+# 0.96
+cd-hit-est \
+-i Merged_Reads/megahit/merged_sequences_GRCH38_p0.fasta \
+-o Cluster_CDHIT/megahit/clustered_GRCH38_p0_n5_i0.96.fasta \
+-c 0.96 \
+-n 5 \
+-T $SLURM_NPROCS  
+
+
+# masurca -----
+
+# 0.9
+cd-hit-est \
+-i Merged_Reads/masurca/merged_sequences_GRCH38_p0.fasta \
+-o Cluster_CDHIT/masurca/clustered_GRCH38_p0_n5_i0.9.fasta \
+-c 0.9 \
+-n 5 \
+-T $SLURM_NPROCS 
+
+# 0.92
+cd-hit-est \
+-i Merged_Reads/masurca/merged_sequences_GRCH38_p0.fasta \
+-o Cluster_CDHIT/masurca/clustered_GRCH38_p0_n5_i0.92.fasta \
+-c 0.92 \
+-n 5 \
+-T $SLURM_NPROCS 
+
+# 0.94
+cd-hit-est \
+-i Merged_Reads/masurca/merged_sequences_GRCH38_p0.fasta \
+-o Cluster_CDHIT/masurca/clustered_GRCH38_p0_n5_i0.94.fasta \
+-c 0.94 \
+-n 5 \
+-T $SLURM_NPROCS 
+
+# 0.96
+cd-hit-est \
+-i Merged_Reads/masurca/merged_sequences_GRCH38_p0.fasta \
+-o Cluster_CDHIT/masurca/clustered_GRCH38_p0_n5_i0.96.fasta \
+-c 0.96 \
+-n 5 \
+-T $SLURM_NPROCS 
+
+
+# CHM13 -----
+
+# megahit
+
+# 0.9
+cd-hit-est \
+-i Merged_Reads/megahit/merged_sequences_CHM13.fasta \
+-o Cluster_CDHIT/megahit/clustered_CHM13_n5_i0.9.fasta \
+-c 0.9 \
+-n 5 \
+-T $SLURM_NPROCS 
+
+# 0.92
+cd-hit-est \
+-i Merged_Reads/megahit/merged_sequences_CHM13.fasta \
+-o Cluster_CDHIT/megahit/clustered_CHM13_n5_i0.92.fasta \
+-c 0.92 \
+-n 5 \
+-T $SLURM_NPROCS 
+
+# 0.94
+cd-hit-est \
+-i Merged_Reads/megahit/merged_sequences_CHM13.fasta \
+-o Cluster_CDHIT/megahit/clustered_CHM13_n5_i0.94.fasta \
+-c 0.94 \
+-n 5 \
+-T $SLURM_NPROCS 
+
+# 0.96
+cd-hit-est \
+-i Merged_Reads/megahit/merged_sequences_CHM13.fasta \
+-o Cluster_CDHIT/megahit/clustered_CHM13_n5_i0.96.fasta \
+-c 0.96 \
+-n 5 \
+-T $SLURM_NPROCS 
+
+
+# masurca -----
+
+# 0.9
+cd-hit-est \
+-i Merged_Reads/masurca/merged_sequences_CHM13.fasta \
+-o Cluster_CDHIT/masurca/clustered_CHM13_n5_i0.9.fasta \
+-c 0.9 \
+-n 5 \
+-T $SLURM_NPROCS 
+
+# 0.92
+cd-hit-est \
+-i Merged_Reads/masurca/merged_sequences_CHM13.fasta \
+-o Cluster_CDHIT/masurca/clustered_CHM13_n5_i0.92.fasta \
+-c 0.92 \
+-n 5 \
+-T $SLURM_NPROCS 
+
+# 0.94
+cd-hit-est \
+-i Merged_Reads/masurca/merged_sequences_CHM13.fasta \
+-o Cluster_CDHIT/masurca/clustered_CHM13_n5_i0.94.fasta \
+-c 0.94 \
+-n 5 \
+-T $SLURM_NPROCS 
+
+# 0.96
+cd-hit-est \
+-i Merged_Reads/masurca/merged_sequences_CHM13.fasta \
+-o Cluster_CDHIT/masurca/clustered_CHM13_n5_i0.96.fasta \
+-c 0.96 \
+-n 5 \
+-T $SLURM_NPROCS 

annotation.nf

@@ -19,8 +19,8 @@ params.outputDir              = "./results"    /*output folder, must specify pat
 params.assembler              = 'megahit'      /*options: megahit|masurca. Default is masurca*/
 
 /*Parameters for cdhit */
-params.cdhit_identity         = '0.9'          /*proportion of idenitity for clustering using cdhit. Default is 0.9*/
-params.cdhit_wordsize         = '7'            /*word size for cdhit. Default is 7*/
+params.cdhit_identity         = '0.99'          /*proportion of idenitity for clustering using cdhit. Default is 0.9*/
+params.cdhit_wordsize         = '5'            /*word size for cdhit. Default is 5*/
 
 /*Stage*/
 stage = "annotation"

filter.nf

@@ -37,7 +37,7 @@ params.blastr_filter_qcov     = '95'           /*filtering cut off for query cov
 
 /*Parameters for cdhit */
 params.cdhit_identity         = '0.9'          /*proportion of idenitity for clustering. Default is 0.9*/
-params.cdhit_wordsize         = '7'            /*word size for cdhit. Default is 7*/
+params.cdhit_wordsize         = '5'            /*word size for cdhit. Default is 5*/
 
 /*Stage*/
 stage = "filter"

other/test_seqkit.sh

@@ -37,14 +37,48 @@ assembly="megahit"
 module purge
 module load BLAST+/2.10.1-IGB-gcc-8.2.0
 
-# blast to compare two sequences --- 
-for j in pipeline_testing/merged_AFR/${assembly}/*seqkit_kn*.fasta
+# # blast to compare two sequences --- 
+# for j in pipeline_testing/merged_AFR/${assembly}/*seqkit_kn*.fasta
+# do 
+# start=`date +%s` # capture start time
+# newj=$(echo $(basename ${j}) | sed 's/.fasta//g')
+# echo "BLAST to reference genome is starting for ${j}"
+# blastn -db ../GRCh38/GRCh38_full_analysis_set_plus_decoy_hla.fa \
+# -query ${j} \
+#  -outfmt "6 qseqid sseqid stitle pident \
+#     length evalue qcovs bitscore mismatch \
+#     gapopen qstart qend qlen sstart send slen" \
+# -perc_identity 100 \
+# -max_target_seqs 5 \
+# -max_hsps 10 \
+# -evalue 1e-5 \
+# -num_threads $SLURM_NPROCS \
+# -out pipeline_testing/BLAST/${assembly}/${newj}_blast_ref.0.tsv
+
+# # Create headers for the blast file -----
+# echo -e "qseqid,sseqid,stitle,pident,length,evalue,qcovs,bitscore,mismatch,gapopen,qstart,qend, \
+# qlen,sstart,send,slen" | tr ',' '\t'| cat - pipeline_testing/BLAST/${assembly}/${newj}_blast_ref.0.tsv \
+# > pipeline_testing/BLAST/${assembly}/${newj}_blast_ref.tsv
+
+# rm pipeline_testing/BLAST/${assembly}/${newj}_blast_ref.0.tsv
+
+# end=`date +%s`
+# runtime=$((end-start))
+# runtime=$( echo "scale=2;$((end-start)) / 60" | bc )
+# echo "It took $runtime minutes to run blast ref genome for ${j}"
+
+# done
+
+
+# blast to compare two sequences ---
+for i in pipeline_testing/merged_AFR/${assembly}/3.seqkit_kn_filtered_All_AFR_megahit.fasta
 do 
-start=`date +%s` # capture start time
-newj=$(echo $(basename ${j}) | sed 's/.fasta//g')
-echo "BLAST to reference genome is starting for ${j}"
-blastn -db ../GRCh38/GRCh38_full_analysis_set_plus_decoy_hla.fa \
--query ${j} \
+start=`date +%s` # capture start time 
+newi=$(echo $(basename ${i}) | sed 's/.fasta//g')
+echo "Subject-Query BLAST is starting for ${i}"
+blastn \
+-query ${i} \
+-subject /home/groups/h3abionet/RefGraph/results/2022-01-30-HGSVCv2-truthset/data/insertions.fasta \
  -outfmt "6 qseqid sseqid stitle pident \
     length evalue qcovs bitscore mismatch \
     gapopen qstart qend qlen sstart send slen" \
@@ -53,19 +87,19 @@ blastn -db ../GRCh38/GRCh38_full_analysis_set_plus_decoy_hla.fa \
 -max_hsps 10 \
 -evalue 1e-5 \
 -num_threads $SLURM_NPROCS \
--out pipeline_testing/BLAST/${assembly}/${newj}_blast_ref.0.tsv
+-out pipeline_testing/BLAST/${assembly}/${newi}_blast_compare.0.tsv
 
 # Create headers for the blast file -----
 echo -e "qseqid,sseqid,stitle,pident,length,evalue,qcovs,bitscore,mismatch,gapopen,qstart,qend, \
-qlen,sstart,send,slen" | tr ',' '\t'| cat - pipeline_testing/BLAST/${assembly}/${newj}_blast_ref.0.tsv \
-> pipeline_testing/BLAST/${assembly}/${newj}_blast_ref.tsv
+qlen,sstart,send,slen" | tr ',' '\t'| cat - pipeline_testing/BLAST/${assembly}/${newi}_blast_compare.0.tsv \
+> pipeline_testing/BLAST/${assembly}/${newi}_blast_compare.tsv
 
-rm pipeline_testing/BLAST/${assembly}/${newj}_blast_ref.0.tsv
+rm pipeline_testing/BLAST/${assembly}/${newi}_blast_compare.0.tsv
 
 end=`date +%s`
 runtime=$((end-start))
 runtime=$( echo "scale=2;$((end-start)) / 60" | bc )
-echo "It took $runtime minutes to run blast ref genome for ${j}"
+echo "It took $runtime minutes to run blast subject-query for ${i}"
 
 done
 

test_compare_insertions_AFRs.sh

@@ -84,7 +84,7 @@ cat filter/Final-Filtered/${assembly}/*_CHM13_filter.final.fasta \
 
 # CD-HIT in annotation by blast GRCH38.p0 -----
 rm pipeline_testing/merged_AFR/${assembly}/7.clustered_*.fasta
-for file in annotation/Cluster_CDHIT/masurca/clustered_*.fasta 
+for file in annotation/Cluster_CDHIT/${assembly}/clustered_*.fasta 
 do
 newfile=$(echo $(basename ${file}) | sed 's/.fasta//g')
 cp -a ${file} pipeline_testing/merged_AFR/${assembly}/7.${newfile}_${assembly}.fasta
@@ -261,7 +261,7 @@ head -n 1 pipeline_testing/SeqKit/${assembly}/All_seqkit_${assembly}.tsv \
 
 # Create headers -----
 echo -e "seqid,format,type,num_seqs,sum_len,min_len,avg_len,max_len,Q1,Q2,Q3,sum_gap,N50,Q20(%),Q30(%)" \
-| tr ',' '\t'| cat - pipeline_testing/Count/${assembly}/shared_sequences_qstats0.tsv \
+| tr ',' '\t' | cat - pipeline_testing/Count/${assembly}/shared_sequences_qstats0.tsv \
 > pipeline_testing/Count/${assembly}/shared_sequences_qstats_${assembly}.tsv
 
 # Count number of pident of 100 and qcov 100 for each file and put into a table ------
@@ -360,7 +360,7 @@ main ()
   echo ""
 	echo "*** RUNNING ${runtype} BLAST COMPARE PIPELINE ***"
 
-quast
+#quast
 #blast_test
 #blast_ref
 count