Mercurial > repos > namhsuya > te_finder

comparison TEfinder.sh @ 0:b81a83c743d3 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author namhsuya
date Tue, 09 Aug 2022 06:58:49 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:b81a83c743d3
1 #!/bin/bash
2 ##
3 ##
4 ## Authors: Vista Sohrab & Dilay Hazal Ayhan
5 ## Date: January 15, 2021
6 ## Description: TEfinder uses discordant reads to detect novel transposable element insertion events in short read paired-end sample sequencing data.
7 ## Software dependencies include bedtools 2.28.0 or later, samtools 1.3 or later, picard 2.0.1 or later
8 ## Required inputs include sample alignment file (.bam|.sam), reference genome FASTA (.fa), reference TE annotation in GFF/GTF or GFF3 (.gff|.gtf), and TEs of interest (.txt)
9 ##
10 ## University of Massachusetts Amherst
11 ##
12 ##
13 ##
14 ##
15
16 set -e
17
18 margs=4
19
20 # Functions
21 function example {
22 echo -e "example: TEfinder -alignment sample.bam -fa reference.fa -gtf TEs.gtf -te List_of_TEs.txt"
23 }
24
25 function help {
26 echo -e "REQUIRED:"
27 echo -e " -alignment, --alignmentFile STR sample reads aligned to reference genome (BAM/SAM file)"
28 echo -e " -fa, --FastaFile STR reference genome FASTA index (FA file)"
29 echo -e " -gtf, --TransposonsInGenome STR reference genome TE annotation (GFF2/GTF file)"
30 echo -e " -te, --TransposonsToSearch STR TE names (single column text file)\n"
31 echo -e "OPTIONAL:"
32 echo -e " -bamo, --DiscordantReads STR BAM output\n"
33 echo -e " -bedo, --bTEinsertions STR TEinsertions BED output\n"
34 echo -e " -gtfo, --gTEinsertions STR TEinsertions GTF output\n"
35 echo -e " -fis, --FragmentInsertSize INT short-read sequencing fragment insert size [400]"
36 echo -e " -picard, --pathToPicardjar STR path to picard tools .jar file [picard.jar]"
37 echo -e " -md, --MaxDistanceForMerge INT maximum distance between reads for bedtools merge [150]"
38 echo -e " -k, --MaxTSDLength INT maximum TE target site duplication (TSD) length [20]"
39 echo -e " -maxHeapMem, --MaxHeapMemory INT java maximum heap memory allocation for picard in Mb [2000]"
40 echo -e " -workingdir, --WorkingDirectory STR working directory name [TEfinder_<Date>]"
41 echo -e " -out, --OutputFormat STR output format as GTF [BED]"
42 echo -e " -outname, --OutputName STR output name prefix added to file names [null]"
43 echo -e " -threads, --Threads INT number of threads for samtools multi-threading [1]"
44 echo -e " -intermed --IntermediateFiles STR keep intermediate files created by pipeline [no]"
45 echo -e " -h, --help prints help\n"
46 example
47 }
48
49 # check if mandatory args are empty
50 function margs_check {
51 if [ $# -lt $margs ]; then
52 echo -e "One or more required parameters are missing."
53 example
54 exit 1 # error
55 fi
56 }
57
58 # main workflow
59 #### : comment out
60 function pipeline() {
61 mkdir ${workingdir}/${line}
62 currdir=${workingdir}/${line}
63 echo -e $(date) " Transposon analysis for "${line}" has started\n"
64
65 grep -P '[^(\w|\d|\-|\_|\#|\.)]'${line}'[^(\w|\d|\-|\_|\#|\.)]' $gtf > ${currdir}/${line}_TE.gff
66 echo -e $(date) " Individual TE GFF has been created for "${line}"\n" ####
67
68 bedtools intersect -abam ${workingdir}/${outname}Alignments.bam -b ${currdir}/${line}_TE.gff -wa > ${currdir}/${line}_MappedReadsToTE.bam
69 echo -e $(date) " Mapped reads to TE via bedtools intersect has been completed for "${line}"\n" ####
70 samtools view -@ $threads ${currdir}/${line}_MappedReadsToTE.bam | \
71 awk -v Ins=`expr $fis \* 10` '{if (($7 != "=") || ($9 > Ins) || ($9 < -Ins)) print $1}' > ${currdir}/${line}_ReadID.txt
72 echo -e $(date) " Identifying discordant read IDs has been completed for "${line}"\n" ####
73
74 # if discordant readID file exists, then continue with remainder of TE analysis
75 if [[ -s ${currdir}/${line}_ReadID.txt ]]
76 then
77 #java $maxHeapMem -jar $picard FilterSamReads I=${workingdir}/${outname}Alignments.bam O=${currdir}/${line}_DiscordantPairs.bam READ_LIST_FILE=${currdir}/${line}_ReadID.txt FILTER=includeReadList WRITE_READS_FILES=false
78 $picard $maxHeapMem FilterSamReads -I ${workingdir}/${outname}Alignments.bam -O ${currdir}/${line}_DiscordantPairs.bam -READ_LIST_FILE ${currdir}/${line}_ReadID.txt -FILTER includeReadList -WRITE_READS_FILES false &>/dev/null
79 echo -e $(date) " Filtering original alignment based on discordant reads IDs is complete for "${line}"\n" ####
80
81 bedtools merge -d $md -S + -c 1 -o count -i ${currdir}/${line}_DiscordantPairs.bam | \
82 awk '{if ($4 > 3) print $0}' > ${currdir}/${line}_plusCluster.bed
83 echo -e $(date) " Primary reads from the + strand have been merged if read count greater than 3 for "${line}"\n" ####
84
85 bedtools merge -d $md -S - -c 1 -o count -i ${currdir}/${line}_DiscordantPairs.bam | \
86 awk '{if ($4 > 3) print $0}' > ${currdir}/${line}_minusCluster.bed
87 echo -e $(date) " Primary reads from the - strand have been merged if read count greater than 3 for "${line}"\n" ####
88
89 # filtering edges piped into bedtools merge (keeping read counts greater than 3 in the line above)
90 ## find the closest minus strand to the plus strand in the cluster
91 ## filter by the distance between the plus and minus clusters - only retain pairs if reads are 0-100 bases away
92 ## if plus strand start is less than minus strand start and plus strand end is less than minus strand end then in proper orientation
93 bedtools closest -d -g ${workingdir}/reference.fa.fai -t first -a ${currdir}/${line}_plusCluster.bed -b ${currdir}/${line}_minusCluster.bed | \
94 awk -v TSD=$k '{if ($9 <= TSD && $9 >= 0) print $0}' | \
95 awk '{if ($2 < $6 && $3 < $7) print $0}' > ${currdir}/${line}_plusminus.bed
96 echo -e $(date) " Filtration of clusters in proper orientation using bedtools closest has been completed for "${line}"\n" ####
97
98 # if plus strand end is greater than minus strand start, then report the pair
99 awk '{if ($3 > $6) print $1"\t"$6"\t"$3"\t"$0}' ${currdir}/${line}_plusminus.bed > ${currdir}/${line}_plusminus_1.bed
100 echo -e $(date) " Overlapping reads TE insertions reported for "${line}"\n" ####
101
102 #if plus strand end is less than or equal to minus strand start and the region in between is less than a user-defined value k, report the pair
103 awk -v TSD=$k '{if ($3 <= $6 && $6 - $3 < TSD) print $1"\t"$3 - 1"\t"$6 + 1"\t"$0}' ${currdir}/${line}_plusminus.bed > \
104 ${currdir}/${line}_plusminus_2.bed
105 echo -e $(date) " Non-overlapping reads TE insertions reported for "${line}"\n" ####
106
107 #combine reported TE insertions
108 cat ${currdir}/${line}_plusminus_1.bed ${currdir}/${line}_plusminus_2.bed | \
109 awk -v TEname=$line '{$0=TEname"\t"$0}1' | sort -k 1 | sort -k 2 > \
110 ${currdir}/${line}_insertionRegion.txt
111
112
113 cat ${currdir}/${line}_insertionRegion.txt >> ${workingdir}/insertions.txt
114 echo -e $(date) " TE insertions for "${line}" have been reported.\n" ####
115
116
117 echo -e $(date) " Transposon named "${line}" is processed.\n"
118 else
119 echo -e $(date) " Transposon named "${line}" is processed. No discordant reads found.\n"
120 rm -r ${currdir}
121 fi
122 }
123
124 # functions end
125
126 # get arguments
127 fa=
128 alignment=
129 gtf=
130 te=
131 out=
132 intermed=
133 bamo="DiscordantReads.bam"
134 bedo="TEinsertions.bed"
135 gtfo="TEinsertions.gtf"
136 outname=""
137 fis=400
138 picard="picard"
139 maxHeapMem=-Xmx2000m
140 md=150
141 k=20
142 d=$(date +%Y%m%d%H%M%S)
143 # workingdir=TEfinder_${d}
144 workingdir="TEfinder"
145 threads=1
146
147 while [ "$1" != "" ];
148 do
149 case $1 in
150 -fa | --FastaFile )
151 shift
152 fa=$1 ;;
153 -alignment | --alignmentFile )
154 shift
155 alignment=$1 ;;
156 -gtf | --TransposonsInGenome )
157 shift
158 gtf=$1 ;;
159 -te | --TransposonsToSearch )
160 shift
161 te=$1 ;;
162 -fis | --FragmentInsertSize )
163 shift
164 fis=$1 ;;
165 -picard | --pathToPicardjar )
166 shift
167 picard=$1 ;;
168 -md | --MaxDistanceForMerge )
169 shift
170 md=$1 ;;
171 -k | --MaxTSDLength )
172 shift
173 k=$1 ;;
174 -bamo | --DiscordantReads )
175 shift
176 bamo=$1 ;;
177 -bedo | --bTEinsertions )
178 shift
179 bedo=$1 ;;
180 -gtfo | --gTEinsertions )
181 shift
182 gtfo=$1 ;;
183 -maxHeapMem | --MaxHeapMemory )
184 shift
185 maxHeapMem="-Xmx"$1"m" ;;
186 -workingdir | --WorkingDirectory )
187 shift
188 workingdir=$1 ;;
189 -out | --OutputFormat )
190 shift
191 out=$1 ;;
192 -outname | --OutputName )
193 shift
194 outname=$1 ;;
195 -threads | --Threads )
196 shift
197 threads=$1 ;;
198 -intermed | --IntermediateFiles )
199 shift
200 intermed=$1 ;;
201 -h | --help )
202 help
203 exit;;
204 *) # error
205 echo "TEfinder: illegal option $1"
206 example
207 exit 1 ;;
208 esac
209 shift
210 done
211 margs_check $fa $alignment $gtf $te
212
213 # main
214
215 mkdir ${workingdir}
216
217 # remove empty lines from user provided TE list if present
218 sed '/^$/d' $te > ${workingdir}"/userTE_noEmptyLines.txt"
219
220 # create output file
221 printf "%s\t" "track name=TEfinder" "type=bedDetail" "description=FR:forward read, RR:reverse read, InsRegion:insertion region start and end positions, FILTER:comma separated filters" > ${workingdir}/${outname}TEinsertions.bed
222 printf "\n" >> ${workingdir}/${outname}TEinsertions.bed
223
224 # create fasta index (fai) file
225 cp $fa ${workingdir}/reference.fa
226 samtools faidx ${workingdir}/reference.fa
227
228 # sort alignment input
229 samtools sort -@ $threads -o ${workingdir}/alignmentInput.sorted.bam ${alignment}
230 echo -e $(date) " Alignment file sorted successfully.\n"
231
232 # remove secondary and supplementary alignments from sorted bam
233 samtools view -F 2304 -@ $threads -o ${workingdir}/${outname}Alignments.bam ${workingdir}/alignmentInput.sorted.bam
234 echo -e $(date) " Alignments are filtered - secondary and supplementary alignments have been removed. \n"
235
236 # run pipeline for each TE
237 while IFS="" read -r line || [ -n "$line" ]
238 do
239 pipeline &
240 done < ${workingdir}/userTE_noEmptyLines.txt
241 wait
242 echo -e $(date) " All transposons are processed. Finalizing...\n"
243
244 # combine discordant bam files
245 samtools merge -@ $threads -r ${workingdir}/${outname}DiscordantReads.bam ${workingdir}/*/*_DiscordantPairs.bam
246 echo -e $(date) " BAM Output: Discordant pair alignment file is now available.\n"
247 # Sorting by position
248 samtools sort -@ $threads ${workingdir}/${outname}DiscordantReads.bam | samtools view -h -o ${workingdir}/${outname}DiscordantReads.sam
249 grep -v '^@PG' ${workingdir}/${outname}DiscordantReads.sam > ${workingdir}/${outname}DiscordantReadsNoPG.sam
250 rm ${workingdir}/${outname}DiscordantReads.sam
251 samtools view -hb -x "PG" --no-PG --remove-flags "PG" -O BAM ${workingdir}/${outname}DiscordantReadsNoPG.sam -o ${bamo}
252 rm ${workingdir}/${outname}DiscordantReadsNoPG.sam
253
254 # update output BED file with TEfinder results: organize the starting file
255 awk '{print $2"\t"$3"\t"$4"\t"$1"\t"$8+$12"\t.\tFR="$8";RR="$12";InsRegion="$6"-"$11";FILTER="}' ${workingdir}/insertions.txt > ${workingdir}/TEinsertions_putative.bed
256 # find the entries in repeat regions for filtering
257 bedtools intersect -wa -u -a ${workingdir}/TEinsertions_putative.bed -b $gtf > ${workingdir}/TEinsertions_putative_inrepeat.bed
258 # filtering process
259 while IFS="" read -r line || [ -n "$line" ]
260 do
261 #located in repeat region
262 if (grep -Fxq "$line" "${workingdir}/TEinsertions_putative_inrepeat.bed")
263 then
264 line=$line"in_repeat,"
265 fi
266
267 #weak evidence
268 readc=$(echo $line | awk '{print $5}')
269 if (( $readc < 10 ))
270 then
271 line=$line"weak_evidence,"
272 fi
273
274 #strand-biased
275 FR=$(echo $line | grep -o 'FR=[[:digit:]]*' | cut -f2 -d'=')
276 RR=$(echo $line | grep -o 'RR=[[:digit:]]*' | cut -f2 -d'=')
277 var1=$(echo 'e(l('$FR')*1.25)' | bc -l)
278 var2=$(echo 'e(l('$FR')*0.8)' | bc -l)
279
280 if [ $(echo "$RR > $var1" | bc) -eq 1 ] || [ $(echo "$RR < $var2" | bc) -eq 1 ]
281 then
282 line=$line"strand_bias,"
283 fi
284
285 #pass
286 lastchar=${line: -1}
287 if [ $lastchar == "," ]
288 then
289 line=${line::${#line}-1}
290 else
291 line=$line"PASS"
292 fi
293
294 #write to final output
295 printf "%s\n" "$line" >> ${workingdir}/${outname}TEinsertions.bed
296
297 done < ${workingdir}/TEinsertions_putative.bed
298 wait
299 echo -e $(date) " BED Output: TEfinder output BED file is now available.\n"
300 # Sorting
301 # cp ${workingdir}/${outname}TEinsertions.bed ${outo}
302 bedtools sort -chrThenSizeA -i ${workingdir}/${outname}TEinsertions.bed > ${bedo}
303 # cat ${bedo}
304
305 # gtf option - create output GTF files with TEfinder results
306 if [ ! -z "$out" ]
307 then
308 awk 'FNR > 1 {print $1"\tTEfinder\tTIP\t"$2 + 1"\t"$3"\t"$5"\t.\t.\tte_name \""$4"\"; tags \""$7"\""}' ${workingdir}/${outname}TEinsertions.bed > ${workingdir}/${outname}TEinsertions.gtf
309 bedtools sort -chrThenSizeA -i ${workingdir}/${outname}TEinsertions.gtf > ${gtfo}
310 # awk 'FNR > 1 {print $1"\tTEfinder\tTIP\t"$2 + 1"\t"$3"\t"$5"\t.\t.\tte_name \""$4"\"; tags \""$7"\""}' ${bedo} > ${gtfo}
311 # Sorting
312 # cp ${workingdir}/${outname}TEinsertions.gtf ${gtfo}
313 echo -e "\n\n"
314 # cat ${gtfo}
315 # bedtools sort -chrThenSizeA -i ${workingdir}/${outname}TEinsertions.gtf > ${gtfo}
316 echo -e $(date) " GTF Output: TEfinder output GTF file is now available.\n"
317 fi
318
319 # clean working directory
320 if [ -z "$intermed" ]
321 then
322 rm ${workingdir}/TEinsertions_putative.bed ${workingdir}/TEinsertions_putative_inrepeat.bed ${workingdir}/reference.fa ${workingdir}/reference.fa.fai \
323 ${workingdir}/alignmentInput.sorted.bam ${workingdir}/insertions.txt ${workingdir}/${outname}Alignments.bam ${workingdir}/userTE_noEmptyLines.txt
324 rm -r ${workingdir}/*/
325 fi
326
327 if [ `wc -l <${workingdir}/${outname}TEinsertions.bed` -le "1" ]
328 then
329 echo -e $(date) " Error: TEfinder run unsuccessful."
330 else
331 echo -e $(date) " TE insertion output files have been created. TEfinder completed successfully."
332 fi