Skip to content

added humann2 script #1

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Marysteph wants to merge 1 commit into
base: master
Choose a base branch
from
Open

added humann2 script #1

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
52 changes: 52 additions & 0 deletions scripts/HUMAnN2_pipeline.sh
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,52 @@
#!/usr/bin/env bash
#$ -o ./log/
#$ -e ./log/


#STEPS
# run the humann2
# renormalize to copies per million and relative abundance
# map to KO and EC terms
# remove the intermediate files

#specify paths

@tkosciolek tkosciolek Mar 18, 2020

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should the user be chaning anything else besides those paths below (lines 14-17)? if not, please mark this section clearly, so that there is no doubt the user should only change those paths. If there is anything else to be changed, I am a big fan of grouping those things together, so that the user does not need to go through the whole script to change variables.

@Marysteph Marysteph Mar 19, 2020

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

here user only needs to change the paths only. Will automate the workflow and make some of the things clear.


WORKING_DIR=/klaster/scratch/mmaranga/functional_analysis/diabimmune_data/
READ_DIR=${WORKING_DIR}/combined_data
OUTPUT_DIR=${WORKING_DIR}/uniref90_diabimmune_humman2_analysis
MAPPING_FILES=/home/mmaranga/utility_mapping/

cd $WORKING_DIR

# running humann2

@tkosciolek tkosciolek Mar 18, 2020

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

how long does the pipeline take? can I run it on a single computer or do I need to run it on a cluser? If it needs to be run on a cluster, please add a docstring at the top specifying that and include sample parameters for a specific queuing system. If the script can be run on a single computer, but it takes more than a few minutes, I'd consider changing the comments (like #running humann2) into print statements (echo "running humann2...), so that the user knows the script is executing.

@Marysteph Marysteph Mar 19, 2020

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay will make changes on that. Diamond translated search part is very computational intensive. Takes like a day for one sample. Would't be advisable to run the pipeline on a personal computer. Especially if one has multiple files


for file in $READ_DIR/*fastq.gz
do
humann2 -i ${file} -o $OUTPUT_DIR
done

# Merge gene family and abundance files

humann2_join_tables -i $OUTPUT_DIR -o $OUTPUT_DIR/merged_genefamilies.tsv --file_name genefamilies

humann2_join_tables -i $OUTPUT_DIR -o $OUTPUT_DIR/merged_pathabundance.tsv --file_name pathabundance

humann2_join_tables -i $OUTPUT_DIR -o $OUTPUT_DIR/merged_pathcoverage.tsv --file_name pathcoverage

# normalization
#normalizing to copies per million
humann2_renorm_table -i $OUTPUT_DIR/merged_genefamilies.tsv -o $OUTPUT_DIR/normalized_genefamilies_cpm.tsv --units cpm
humann2_renorm_table -i $OUTPUT_DIR/merged_pathabundance.tsv -o $OUTPUT_DIR/normalized_pathabundance_cpm.tsv --units cpm

# normalizing to relative abundance
humann2_renorm_table -i $OUTPUT_DIR/merged_genefamilies.tsv -o $OUTPUT_DIR/normalized_genefamilies_relab.tsv --units relab
humann2_renorm_table -i $OUTPUT_DIR/merged_pathabundance.tsv -o $OUTPUT_DIR/normalized_pathabundance_relab.tsv --units relab

# map to KO
humann2_regroup_table -i $OUTPUT_DIR/normalized_genefamilies_cpm.tsv -o $OUTPUT_DIR/normalized_genefamilies_cpm_KO.tsv -c $MAPPING_FILES/map_ko_uniref90.txt.gz

# map genefamilies abundance to Enzyme Commission (EC) terms
humann2_regroup_table -i $OUTPUT_DIR/normalized_genefamilies_cpm.tsv -o $OUTPUT_DIR/normalized_genefamilies_cpm_EC.tsv -c $MAPPING_FILES/map_level4ec_uniref90.txt.gz

# remove the intermediate files
rm -rf $OUTPUT_DIR/*_humann2_temp

@tkosciolek tkosciolek Mar 18, 2020

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

are those temp files potentially useful for debugging? if yes, then consider not removing them by default. or using a flag to determine if temp files should be kept or removed. For example using the getopts function (e.g. https://stackoverflow.com/questions/14447406/bash-shell-script-check-for-a-flag-and-grab-its-value)

@Marysteph Marysteph Mar 19, 2020

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

temp files can be useful for debugging. They were so huge, especially sam files from bowtie step. They were taking so much space on the cluster. I will consider be keeping them