|
0
|
1 #!/bin/bash
|
|
|
2 set -e; mkdir -p vcf; mkdir -p bam; mkdir -p references; mkdir -p scripts;
|
|
|
3
|
|
|
4 email_address_of_user=$1;
|
|
|
5 name_of_user_computer=$2;
|
|
|
6
|
|
|
7 #read config variables
|
|
|
8 echo "Reading config...." >&2
|
|
|
9 #echo `ls -l`
|
|
|
10 computer_name=`hostname`;
|
|
|
11 if [ $computer_name == $name_of_user_computer ]; then
|
|
|
12 #local_computer
|
|
|
13 source config_devel;
|
|
|
14 else
|
|
|
15 #metacentrum machine
|
|
|
16 source config_meta;
|
|
|
17 fi
|
|
|
18
|
|
|
19 #initally we believe that all files are present
|
|
|
20 ok=true;
|
|
|
21
|
|
|
22 #we get to know if files are provided and if they are in .fasta or .fastq format
|
|
|
23 #first, rename files to shorter versions (fastq->fq, fasta->fa)
|
|
|
24 rename 's/fasta$/fa/' input/*.fasta 2>/dev/null
|
|
|
25 rename 's/fastq$/fq/' input/*.fastq 2>/dev/null
|
|
|
26
|
|
|
27 check_if_exists ()
|
|
|
28 {
|
|
|
29 for ARG in "$@";
|
|
|
30 do
|
|
|
31 if [[ -e "input/${ARG}_1.fq" ]]; then
|
|
|
32 eval ${ARG}_ext=fq
|
|
|
33 elif [[ -e "input/${ARG}_1.fa" ]]; then
|
|
|
34 eval ${ARG}_ext=fa
|
|
|
35 else
|
|
|
36 ok=false;
|
|
|
37 fi
|
|
|
38 done;
|
|
|
39 }
|
|
|
40
|
|
|
41 check_if_paired_end ()
|
|
|
42 {
|
|
|
43 for ARG in "$@";
|
|
|
44 do
|
|
|
45 extension=$(eval "echo \$${ARG}_ext"); #echo $extension;
|
|
|
46
|
|
|
47 if [[ (-e "input/${ARG}_1.${extension}") && (-e "input/${ARG}_2.${extension}") ]]; then
|
|
|
48 eval ${ARG}_paired=paired;
|
|
|
49 else
|
|
|
50 eval ${ARG}_paired=single;
|
|
|
51 fi
|
|
|
52 done;
|
|
|
53 }
|
|
|
54
|
|
|
55 #check if file with reads exist
|
|
|
56 check_if_exists mother father daughter son;
|
|
|
57 #and determine whether they are single or paired-end reads
|
|
|
58 check_if_paired_end mother father daughter son;
|
|
|
59
|
|
|
60 if [ $ok == true ]; then
|
|
|
61 echo "mother reads names:" $mother_ext $mother_paired;
|
|
|
62 echo "father reads names:" $father_ext $father_paired;
|
|
|
63 echo "daughter reads names:" $daughter_ext $daughter_paired;
|
|
|
64 echo "son reads names:" $son_ext $son_paired;
|
|
|
65
|
|
|
66 if [[ $father_ext =~ (fq) ]] && [[ $son_ext =~ (fq) ]]; then
|
|
|
67 #no conversion needed, both fastq files
|
|
|
68 male_ext=fq;
|
|
|
69 else
|
|
|
70 #convert to fasta (we need to merge father and son reads in order to build male reference contigs)
|
|
|
71 #father OR son is fasta
|
|
|
72 echo "converting..";
|
|
|
73 male_ext=fa;
|
|
|
74
|
|
|
75 if [[ father_1.$father_ext =~ (fq) ]]; then
|
|
|
76 #father is fastq and must be converted to fasta
|
|
|
77 fastq_to_fasta -n -Q33 -i input/father_1.$father_ext -o input/father_1.$male_ext;
|
|
|
78 fastq_to_fasta -n -Q33 -i input/father_2.$father_ext -o input/father_2.$male_ext;
|
|
|
79 father_ext=fa;
|
|
|
80
|
|
|
81 else
|
|
|
82 #son is fastq and must be converted to fasta
|
|
|
83 fastq_to_fasta -n -Q33 -i input/son_1.$father_ext -o input/son_1.$male_ext;
|
|
|
84 fastq_to_fasta -n -Q33 -i input/son_2.$father_ext -o input/son_2.$male_ext;
|
|
|
85 son_ext=fa;
|
|
|
86 fi
|
|
|
87 fi
|
|
|
88 #merging of father and son files
|
|
|
89 echo "concatenating..";
|
|
|
90 cat input/father_1.$father_ext input/son_1.$son_ext >input/father_and_sons_1.$male_ext;
|
|
|
91
|
|
|
92 if [[ $father_paired =~ (paired) ]] && [[ $son_paired =~ (paired) ]]; then
|
|
|
93 cat input/father_2.$father_ext input/son_2.$son_ext >input/father_and_sons_2.$male_ext;
|
|
|
94 fi
|
|
|
95
|
|
|
96
|
|
|
97 #NOW run Trinity to assemble mother contigs
|
|
|
98 #in: mother reads out: reference.fasta
|
|
|
99 PBS_NUM_PPN_used=$(($PBS_NUM_PPN_used - 2));
|
|
|
100 echo "processors used:" $PBS_NUM_PPN_used;
|
|
|
101
|
|
|
102 #<<PREVENT_TRINITY_MOTHER
|
|
|
103 #assemble mother reads in single or paired-end mode
|
|
|
104 if [[ $mother_paired =~ (paired) ]]; then
|
|
|
105 $trinity --CPU $PBS_NUM_PPN_used --seqType $mother_ext --JM 20G --left input/mother_1.$mother_ext --right input/mother_2.$mother_ext --output trinity_out_dir_mother 1>$trinity_stdout 2>&1
|
|
|
106 echo "$trinity --CPU $PBS_NUM_PPN_used --seqType $mother_ext --JM 20G --left input/mother_1.$mother_ext --right input/mother_2.$mother_ext --output trinity_out_dir_mother 1>$trinity_stdout 2>&1";
|
|
|
107 else
|
|
|
108 $trinity --CPU $PBS_NUM_PPN_used --seqType $mother_ext --JM 20G --single input/mother_1.$mother_ext --output trinity_out_dir_mother 1>$trinity_stdout 2>&1
|
|
|
109 echo "$trinity --CPU $PBS_NUM_PPN_used --seqType $mother_ext --JM 20G --single input/mother_1.$mother_ext --output trinity_out_dir_mother 1>$trinity_stdout 2>&1";
|
|
|
110 fi
|
|
|
111
|
|
|
112 #retrieve assembled contigs and name them as reference.fasta
|
|
|
113 mv trinity_out_dir_mother/Trinity.fasta reference_pro.fasta
|
|
|
114 perl -pe 's/(>[\S]+).+/$1/' reference_pro.fasta >references/reference.fasta; rm reference_pro.fasta;
|
|
|
115 #PREVENT_TRINITY_MOTHER
|
|
|
116
|
|
|
117 grep ">" references/reference.fasta | sed 's/>//g' >references/reference_contig_names
|
|
|
118
|
|
|
119 #now call to_get_vcf_files.sh to map reads to reference contigs
|
|
|
120 ./scripts/to_get_vcf_files.sh references/reference.fasta mother_1.$mother_ext mother_2.$mother_ext mother $mother_paired
|
|
|
121 echo "Mother mapping successful.";
|
|
|
122 ./scripts/to_get_vcf_files.sh references/reference.fasta father_1.$mother_ext father_2.$mother_ext father $father_paired
|
|
|
123 echo "Father mapping successful.";
|
|
|
124 ./scripts/to_get_vcf_files.sh references/reference.fasta daughter_1.$mother_ext daughter_2.$mother_ext daughter $daughter_paired
|
|
|
125 echo "Daughter mapping successful.";
|
|
|
126 ./scripts/to_get_vcf_files.sh references/reference.fasta son_1.$mother_ext son_2.$mother_ext son $son_paired
|
|
|
127 echo "Son mapping successful.";
|
|
|
128
|
|
|
129 wait;
|
|
|
130
|
|
|
131 #NOW run Trinity to assemble male contigs
|
|
|
132 #in: male reads out: reference_male.fasta
|
|
|
133 #<<PREVENT_TRINITY_MALES
|
|
|
134 if [[ $father_paired =~ (paired) ]] && [[ $son_paired =~ (paired) ]]; then
|
|
|
135 $trinity --CPU $PBS_NUM_PPN_used --seqType $male_ext --JM 20G --left input/father_and_sons_1.$male_ext --right input/father_and_sons_2.$male_ext --output trinity_out_dir_male 1>$trinity_stdout 2>&1
|
|
|
136 echo "$trinity --CPU $PBS_NUM_PPN_used --seqType $male_ext --JM 20G --left input/father_and_sons_1.$male_ext --right input/father_and_sons_2.$male_ext --output trinity_out_dir_male 1>$trinity_stdout 2>&1";
|
|
|
137 else
|
|
|
138 $trinity --CPU $PBS_NUM_PPN_used --seqType $male_ext --JM 20G --single input/father_and_sons_1.$male_ext --output trinity_out_dir_male 1>$trinity_stdout 2>&1
|
|
|
139 echo "$trinity --CPU $PBS_NUM_PPN_used --seqType $male_ext --JM 20G --single input/father_and_sons_1.$male_ext --output trinity_out_dir_male 1>$trinity_stdout 2>&1";
|
|
|
140 fi
|
|
|
141
|
|
|
142 #retrieve assembled contigs and name them as reference_male.fasta
|
|
|
143 mv trinity_out_dir_male/Trinity.fasta reference_male_pro.fasta
|
|
|
144 perl -pe 's/(>[\S]+).+/$1/' reference_male_pro.fasta >references/reference_male.fasta; rm reference_male_pro.fasta;
|
|
|
145 #PREVENT_TRINITY_MALES
|
|
|
146
|
|
|
147 grep ">" references/reference_male.fasta | sed 's/>//g' >references/male_contig_names
|
|
|
148
|
|
|
149 #now call to_get_vcf_files.sh to map reads to reference contigs
|
|
|
150 ./scripts/to_get_Y_variants.sh references/reference_male.fasta mother_1.$mother_ext mother_2.$mother_ext mother $mother_paired
|
|
|
151 echo "Mother mapping to males successful.";
|
|
|
152 ./scripts/to_get_Y_variants.sh references/reference_male.fasta father_1.$mother_ext father_2.$mother_ext father $father_paired
|
|
|
153 echo "Father mapping to males successful.";
|
|
|
154 ./scripts/to_get_Y_variants.sh references/reference_male.fasta daughter_1.$mother_ext daughter_2.$mother_ext daughter $daughter_paired
|
|
|
155 echo "Daughter mapping to males successful.";
|
|
|
156 ./scripts/to_get_Y_variants.sh references/reference_male.fasta son_1.$mother_ext son_2.$mother_ext son $son_paired
|
|
|
157 echo "Son mapping to males successful.";
|
|
|
158
|
|
|
159 #index bam files
|
|
|
160
|
|
|
161 index_bam_file ()
|
|
|
162 {
|
|
|
163 for ARG in "$@";
|
|
|
164 do
|
|
|
165 if [ -e "bam/aln_cleaned_sorted_deduplicated_${ARG}.bam.bai" ];
|
|
|
166 then
|
|
|
167 echo "Index already exists. Not being created again.";
|
|
|
168 else
|
|
|
169 samtools index bam/aln_cleaned_sorted_deduplicated_${ARG}.bam &
|
|
|
170 fi
|
|
|
171 done;
|
|
|
172 }
|
|
|
173
|
|
|
174 #index bam files
|
|
|
175 index_bam_file mother father daughter son mother_male father_male daughter_male son_male;
|
|
|
176 wait;
|
|
|
177 echo "Indexes successfully created."
|
|
|
178
|
|
|
179 echo -e "Starting analysis";
|
|
|
180
|
|
|
181 #A variant X1X1 X2Y
|
|
|
182 dir="A/";
|
|
|
183 bash -x ./A.sh $email_address_of_user;
|
|
|
184 perl scripts/sorting.pl A/;
|
|
|
185
|
|
|
186 #Y variant
|
|
|
187 dir="Y/";
|
|
|
188 bash -x ./Y.sh $email_address_of_user;
|
|
|
189
|
|
|
190 wait;
|
|
|
191
|
|
|
192 #sending results (not BAM files)
|
|
|
193 mutt -s "LINKYX: Results of computation for Y-linked genes" $email_address_of_user -a Y_results_sequences.txt < messages_to_user/message_results_sentY;
|
|
|
194 mutt -s "LINKYX: Results of computation for X-linked genes" $email_address_of_user -a X_results_sequences.txt < messages_to_user/message_results_sentX;
|
|
|
195
|
|
|
196
|
|
|
197 echo -e "Done."
|
|
|
198 else
|
|
|
199 echo "not all input files are present in proper form";
|
|
|
200 mutt -s "LINKYX: Not all input files are present in proper form." $email_address_of_user < messages_to_user/message_input_files_error
|
|
|
201 fi
|
|
|
202
|
|
|
203
|
|
|
204
|
|
|
205
|
