-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsarscov2seq
137 lines (107 loc) · 3.76 KB
/
sarscov2seq
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
#! /usr/bin/env bash
#
#
####################################################
#Global Settings
####################################################
#####Path to Ref directory
ref_dir="$HOME/SARSCOV2seq/Ref/"
script_dir="$HOME/SARSCOV2seq/Scripts/"
date=$(date '+%Y-%m-%d')
threads="8"
##### ivar consensus thresholds
freq="0.9"
qual="20"
depth="10"
##### 29903 bp coronavirus 2 isolate Wuhan-Hu-1 MN908947.3 18-MAR-2020
reference="coronavirus2_wuhan-hu-1_MN908947.3_2020-04-29"
reference_fasta="${ref_dir}${reference}.fasta"
reference_genes="${ref_dir}${reference}_genes.txt"
##### Settings for MTBseq
### do not do quality recalibration, resistance and category annotation, set to "none"
variant_list="none"
resistance_variants_info="none"
resistance_regions_info="none"
gene_categories_info="none"
##### Amplicon target regions
###ARTIC sars-cov-2 amplicon panel
#ARTIC sars-cov-2 amplicon panel v3
#target_list="${ref_dir}targets_list_ARTIC-amplicon_sars-cov-2_2020-06-18.txt"
#ARTIC sars-cov-2 amplicon panel v4.1
#target_list="${ref_dir}targets_list_ARTIC-amplicon_sars-cov-2_4-1_2022-02-21.txt"
###BEDfiles for primerregions
###artic v3
#primer_bed="${ref_dir}artic_nCoV2019_400bp_primers.bed"
###artic v4.1
primer_bed="${ref_dir}artic_sars-cov-2_4-1_400bp_primers.bed"
###primer pair information for ivar trim -f to remove reads spanning 2 amplicons
###artic v3
#primer_pairs="${ref_dir}artic_nCoV2019_400bp_primers_ivar_informationfile.txt"
###artic v4.1
primer_pairs="${ref_dir}artic_sars-cov-2_v4-1_400bp_primers_ivar_informationfile.txt"
##################################################
# Necessary folder structure
##################################################
mkdir Fastq
mkdir Consensus
mkdir Called_TrueCodon
##################################################
# sars-cov-2 pipeline
##################################################
eval "$(conda shell.bash hook)"
conda activate
conda activate mtbseq
### BWA mapping
MTBseq -step TBbwa --ref $reference --threads $threads
### refine (merged) bam files
MTBseq -step TBrefine --ref $reference --threads $threads
conda deactivate
rmdir Groups
rmdir Joint
rmdir Amend
rmdir Classification
cd GATK_Bam
mv *gatk.bam ../
mv *gatk.bai ../
cd ..
conda activate ivar
###trim primer, but will also trim sliding window 4 below Q20, remove reads shorter 30bp and unmapped.
bash "${script_dir}starter_ivar_trimming.v04.sh" $primer_bed $threads $primer_pairs
### consensus
bash "${script_dir}starter_ivar_consensus.v04.sh" $freq $qual $depth
conda deactivate
### calculate N-content of consensus *.fa
bash "${script_dir}starter_fasta_n_content.sh"
###activate conda environment
conda activate pangolin
##calculate sequence identity *.fa
bash "${script_dir}starter_sequence_identity.sh" $reference_fasta
###get variant type with pangolin
cat *.fa > ""$date"_all.fasta"
pangolin --update
pangolin --outfile ""$date"_lineage_report.tsv" -t $threads ""$date"_all.fasta"
conda deactivate
mv *gatk.bam ./GATK_Bam/
mv *gatk.bai ./GATK_Bam/
conda activate mtbseq
MTBseq -step TBpile --ref $reference --threads $threads
MTBseq -step TBlist --ref $reference --threads $threads
MTBseq -step TBvariants --ref $reference --threads $threads
MTBseq -step TBstats --ref $reference --threads $threads
conda deactivate
conda deactivate
### combines multiline variants into a single line, recalculates the AA exchange if in a gene.
cd Called
mv *_position_variants_*.tab ../
cd ..
perl "${script_dir}"mtbseq_double_mut.v02.pl *_position_variants_*.tab
mv *_position_variants_*.tab ./Called/
##################################################
# Move files
##################################################
mv *.fastq.gz ./Fastq/
mv *variants*true* ./Called_TrueCodon
mv *.fa ./Consensus/
mv *.qual.txt ./Consensus/
###################################################
exit