annotate TEfinder @ 0:b81a83c743d3 draft default tip

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