-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathSingh_et_al_genome_annotation_scripts.txt
406 lines (325 loc) · 22.9 KB
/
Singh_et_al_genome_annotation_scripts.txt
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
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
#Scripts for the annotation of D.chrysippus2.2 genome - Singh et al. ***
#annotation of each danaus assembly:
#1. producing a protein dataset - combining broad arthropod set with D. plexippus
#2. masking the genome using repeatmodeler/masker
#3. carry out annotation with BRAKER2
###########################################################
# 1. produce protein dataset:
###########################################################
#this comprises of two D.plexippus protein sets
#FIRST
#plexippus proteome set from assembly: GCF_009731565.1
#https://www.ncbi.nlm.nih.gov/assembly/GCF_009731565.1
grep ">" /data/martin/genomics/analyses/Danaus_genome/MB181_trio/annotation/protein_set/uniprot-proteome_UP000596680.fasta | wc -l
#17926
#SECOND
#plexippus amino acid protein set from assembly: GCA_018135715.1
#https://www.ncbi.nlm.nih.gov/assembly/GCA_018135715.1/#/def
#annotation at: https://zenodo.org/record/4470132#.YRKWDRNKhhE
#convert extract amino acid sequences into .aa.fasta file
sconda /ceph/users/amackintosh/.conda/envs/gt/
# tidy up the braker.gff3
sed 's/:=//g' dplex_mex.gff3 > intermediate.gff3
sed 's/mstrg/MSTRG/g' intermediate.gff3 > intermediate_upper.gff3
gt gff3 -sort -tidy -retainids -fixregionboundaries intermediate_upper.gff3 > dplex_mex.tidy.gff3
# extract CDS sequences
gt extractfeat -type CDS -join -retainids -seqfile dplex_mex.fa -matchdescstart -o dplex_mex.gt.cds.fasta dplex_mex.tidy.gff3
# extract protein sequences
gt extractfeat -type CDS -translate -join -retainids -seqfile dplex_mex.fa -matchdescstart -o dplex_mex.gt.aa.fasta dplex_mex.tidy.gff3
# get stats
gt stat -genelengthdistri -genescoredistri -exonlengthdistri -exonnumberdistri -intronlengthdistri -cdslengthdistri -o dplex_mex.gt.stats.txt dplex_mex.tidy.gff3
#COMBINE
cat /data/martin/genomics/analyses/Danaus_genome/MB181_trio/annotation/protein_set/uniprot-proteome_UP000596680.fasta ./dplex_mex.gt.aa.fasta > /data/martin/genomics/analyses/Danaus_genome/MB181_trio/annotation/protein_set/dplex2_uniprot-proteome_UP000596680_and_dplex_mex.fasta
grep ">" /data/martin/genomics/analyses/Danaus_genome/MB181_trio/annotation/protein_set/dplex2_uniprot-proteome_UP000596680_and_dplex_mex.fasta | wc -l
#35933
###########################################################
# 2. repeat annotation:
###########################################################
/ceph/software/utilities/sge/qlogin -pe smp64 64 -N gen_an -l h=bigfoot
mkdir /scratch/rdekayne/repeats
cd /scratch/rdekayne/repeats
sconda repeats
#build repeat library from chrysippus assembly
/ceph/software/repeatmodeler/RepeatModeler.v2.0.1/BuildDatabase -name danaus_chrysippus_2.2 -engine ncbi /data/martin/genomics/analyses/Danaus_genome/Dchry2/Dchry2.2.fa
/ceph/software/repeatmodeler/RepeatModeler.v2.0.1/RepeatModeler -database danaus_chrysippus_2.2 -pa 50 -LTRStruct
#get other lepidoptera repeats
perl /ceph/software/repeatmasker/RepeatMasker-4.1.0/util/queryRepeatDatabase.pl -species Lepidoptera | grep -v "Species:" > Lepidoptera.Repbase.repeatmasker
grep ">" Lepidoptera.Repbase.repeatmasker | wc -l
#1253
grep ">" danaus_chrysippus_2.2-families.fa | wc -l
1653
#combine
cat Lepidoptera.Repbase.repeatmasker danaus_chrysippus_2.2-families.fa > Lepidoptera_and_danaus_chrysippus2.2.repeatmasker
grep ">" Lepidoptera_and_danaus_chrysippus2.2.repeatmasker | wc -l
#2906
#carry out softmasking of assembly for BRAKER2
cp /data/martin/genomics/analyses/Danaus_genome/Dchry2/Dchry2.2.fa .
#/ceph/software/repeatmasker/RepeatMasker-4.1.0/RepeatMasker -e rmblast -pa 48 -s -lib Lepidoptera_and_danaus_chrysippus2.2.repeatmasker -dir /scratch/rdekayne/annotate/masked_output -xsmall -gff ./Dchry2.2.fa
#get repeat landscape to view expansion of gene families:
/ceph/software/repeatmasker/RepeatMasker-4.1.0/RepeatMasker -e rmblast -pa 48 -s -a -xsmall -gccalc -lib ./Lepidoptera_and_danaus_chrysippus2.2.repeatmasker ./Dchry2.2.fa
/ceph/software/repeatmasker/RepeatMasker-4.1.0/util/calcDivergenceFromAlign.pl -s test.divsum Dchry2.2.fa.align
/ceph/software/repeatmasker/RepeatMasker-4.1.0/util/createRepeatLandscape.pl -div test.divsum -g 354020300 > test_danaus_output.html
#run again without CpG i.e. using -noCpGMod (addressing reviewer comment)
cp Dchry2.2.fa.align /data/martin/genomics/analyses/Danaus_genome/Dchry2/braker_annotation_RDK_Dchry2_2/repeatlandscape_review
/ceph/software/repeatmasker/RepeatMasker-4.1.0/util/calcDivergenceFromAlign.pl -s test2.divsum -noCpGMod Dchry2.2.fa.align
/ceph/software/repeatmasker/RepeatMasker-4.1.0/util/createRepeatLandscape.pl -div test2.divsum -g 354020300 > test_danaus_output2.html
#use this output .html for plotting
#make hardmasked too:
sed -e '/^>/! s/[[:lower:]]/N/g' Dchry2.2.fa.masked > Dchry2.2.fa.hardmasked
#and now dplex 1
cp /data/martin/genomics/analyses/Danaus_genome/Dplex_V4_HiC/Dapl_Zhan_v3_HiC.RN.fasta .
awk 'BEGIN{FS=" "}{if(!/>/){print toupper($0)}else{print $1}}' Dapl_Zhan_v3_HiC.RN.fasta > Dapl_Zhan_v3_HiC.RN.uppercase.fasta
/ceph/software/repeatmasker/RepeatMasker-4.1.0/RepeatMasker -e rmblast -pa 48 -s -a -xsmall -gccalc -lib ./Lepidoptera_and_danaus_chrysippus2.2.repeatmasker ./Dapl_Zhan_v3_HiC.RN.uppercase.fasta
#and now dplex 2
cp /data/martin/genomics/analyses/Danaus_genome/Dplex_mex/dplex_mex.fa .
awk 'BEGIN{FS=" "}{if(!/>/){print toupper($0)}else{print $1}}' dplex_mex.fa > dplex_mex.uppercase.fa
/ceph/software/repeatmasker/RepeatMasker-4.1.0/RepeatMasker -e rmblast -pa 48 -s -a -xsmall -gccalc -lib ./Lepidoptera_and_danaus_chrysippus2.2.repeatmasker ./dplex_mex.uppercase.fa
###########################################################
# 3. BRAKER2 annotation + tidy:
###########################################################
#and annotate with braker gff output inputting the softmasked assembly and the protein set from the two D.plex assemblies
sconda annotation
/ceph/software/conda/envs/annotation/bin/braker.pl --cores 48 --gff3 --species=danaus_chrysippus.2.2_braker_a001 --workingdir=/scratch/rdekayne/Dchr_annotation --softmasking --genome=/scratch/rdekayne/repeats/Dchry2.2.fa.masked --prot_seq=/data/martin/genomics/analyses/Danaus_genome/MB181_trio/annotation/protein_set/dplex2_uniprot-proteome_UP000596680_and_dplex_mex.fasta --prg=gth --gth2traingenes --trainFromGth
# tidy up the braker.gff3
sconda /ceph/users/amackintosh/.conda/envs/gt/
gt gff3 -sort -tidy -retainids -fixregionboundaries augustus.hints.gff3 > ./Dchry2.2.fa.masked_annotation_a001.sequences.tidy.gff3
# extract CDS sequences
gt extractfeat -type CDS -join -retainids -seqfile /scratch/rdekayne/repeats/Dchry2.2.fa.masked -matchdescstart -o Dchry2.2.fa.masked_annotation_a001.sequences.gt.cds.fasta Dchry2.2.fa.masked_annotation_a001.sequences.tidy.gff3
# extract protein sequences
gt extractfeat -type CDS -translate -join -retainids -seqfile /scratch/rdekayne/repeats/Dchry2.2.fa.masked -matchdescstart -o Dchry2.2.fa.masked_annotation_a001.sequences.gt.aa.fasta Dchry2.2.fa.masked_annotation_a001.sequences.tidy.gff3
# get stats
gt stat -genelengthdistri -genescoredistri -exonlengthdistri -exonnumberdistri -intronlengthdistri -cdslengthdistri -o Dchry2.2.fa.masked_annotation_a001.sequences.gt.stats.txt Dchry2.2.fa.masked_annotation_a001.sequences.tidy.gff3
#reannotate other assemblies
sed 's/\;//g' /scratch/rdekayne/repeats/Dapl_Zhan_v3_HiC.RN.uppercase.fasta.masked > /scratch/rdekayne/repeats/Dapl_Zhan_v3_HiC.RN.uppercase.fasta.renamed.masked
sed -i 's/=//g' /scratch/rdekayne/repeats/Dapl_Zhan_v3_HiC.RN.uppercase.fasta.renamed.masked
/ceph/software/conda/envs/annotation/bin/braker.pl --cores 48 --gff3 --species=danaus_plexv4_braker_a006 --workingdir=/scratch/rdekayne/Dplex_v4_annotation --softmasking --genome=/scratch/rdekayne/repeats/Dapl_Zhan_v3_HiC.RN.uppercase.fasta.renamed.masked --prot_seq=/data/martin/genomics/analyses/Danaus_genome/MB181_trio/annotation/protein_set/dplex2_uniprot-proteome_UP000596680_and_dplex_mex.fasta --prg=gth --gth2traingenes --trainFromGth
gt gff3 -sort -tidy -retainids -fixregionboundaries augustus.hints.gff3 > ./danaus_plexv4_braker_a006.sequences.tidy.gff3
# extract CDS sequences
gt extractfeat -type CDS -join -retainids -seqfile /scratch/rdekayne/repeats/Dapl_Zhan_v3_HiC.RN.uppercase.fasta.renamed.masked -matchdescstart -o danaus_plexv4_braker_a006.sequences.gt.cds.fasta danaus_plexv4_braker_a006.sequences.tidy.gff3
# extract protein sequences
gt extractfeat -type CDS -translate -join -retainids -seqfile /scratch/rdekayne/repeats/Dapl_Zhan_v3_HiC.RN.uppercase.fasta.renamed.masked -matchdescstart -o danaus_plexv4_braker_a006.sequences.gt.aa.fasta danaus_plexv4_braker_a006.sequences.tidy.gff3
# get stats
gt stat -genelengthdistri -genescoredistri -exonlengthdistri -exonnumberdistri -intronlengthdistri -cdslengthdistri -o danaus_plexv4_braker_a006.sequences.gt.stats.txt danaus_plexv4_braker_a006.sequences.tidy.gff3
/ceph/software/conda/envs/annotation/bin/braker.pl --cores 48 --gff3 --species=danaus_plex_mex_braker_a002 --workingdir=/scratch/rdekayne/Dplex_mex_annotation --softmasking --genome=/data/martin/genomics/analyses/Danaus_genome/Dplex_mex/dplex_mex.uppercase.fa.masked --prot_seq=/data/martin/genomics/analyses/Danaus_genome/MB181_trio/annotation/protein_set/dplex2_uniprot-proteome_UP000596680_and_dplex_mex.fasta --prg=gth --gth2traingenes --trainFromGth
# tidy up the braker.gff3
gt gff3 -sort -tidy -retainids -fixregionboundaries augustus.hints.gff3 > ./danaus_plex_mex_braker_a002.sequences.tidy.gff3
# extract CDS sequences
gt extractfeat -type CDS -join -retainids -seqfile /data/martin/genomics/analyses/Danaus_genome/Dplex_mex/dplex_mex.uppercase.fa.masked -matchdescstart -o danaus_plex_mex_braker_a002.sequences.gt.cds.fasta danaus_plex_mex_braker_a002.sequences.tidy.gff3
# extract protein sequences
gt extractfeat -type CDS -translate -join -retainids -seqfile /data/martin/genomics/analyses/Danaus_genome/Dplex_mex/dplex_mex.uppercase.fa.masked -matchdescstart -o danaus_plex_mex_braker_a002.sequences.gt.aa.fasta danaus_plex_mex_braker_a002.sequences.tidy.gff3
# get stats
gt stat -genelengthdistri -genescoredistri -exonlengthdistri -exonnumberdistri -intronlengthdistri -cdslengthdistri -o danaus_plex_mex_braker_a002.sequences.gt.stats.txt danaus_plex_mex_braker_a002.sequences.tidy.gff3
###########################################################
# 4. genome annotation processing:
###########################################################
#get original assemblies
#tidy dplex_v3/4 assembly:
rsync -avzP ./Zhan_v3_HiC.deeptools.recoded.gtf qmaster:/data/martin/genomics/analyses/Danaus_genome/MB181_trio/annotation/annotation_verification/
# tidy up the braker.gff3
gt gtf_to_gff3 Zhan_v3_HiC.deeptools.recoded.gtf > Dplex.gff3
gt gff3 -sort -tidy -retainids -fixregionboundaries Dplex.gff3 > ./Dplex.tidy.gff3
# extract CDS sequences
gt extractfeat -type CDS -join -retainids -seqfile /data/martin/genomics/analyses/Danaus_genome/Dplex_V4_HiC/Dapl_Zhan_v3_HiC.RN.fasta -matchdescstart -o Dplex.tidy.gff3.cds.fasta Dplex.tidy.gff3
# extract protein sequences
gt extractfeat -type CDS -translate -join -retainids -seqfile /data/martin/genomics/analyses/Danaus_genome/Dplex_V4_HiC/Dapl_Zhan_v3_HiC.RN.fasta -matchdescstart -o Dplex.tidy.gff3.gt.aa.fasta Dplex.tidy.gff3
# get stats
gt stat -genelengthdistri -genescoredistri -exonlengthdistri -exonnumberdistri -intronlengthdistri -cdslengthdistri -o Dplex.tidy.gff3.gt.stats.txt Dplex.tidy.gff3
#add introns to 'original' - 'O' assemblies
#dplex1
sed 's/CDS/exon/g' ./Dplex.tidy.gff3 > Dplex.introninput.tidy.gff3
gt gff3 --addintrons Dplex.introninput.tidy.gff3 > Dplex.introns.tidy.gff3
gt stat -intronlengthdistri -o Dplex.introns.stats.txt Dplex.introns.tidy.gff3
#dplex2
gt gff3 --addintrons dplex_mex.tidy.gff3 > dplex_mex.introns.tidy.gff3
gt stat -intronlengthdistri -o dplex_mex.introns.stats.txt dplex_mex.introns.tidy.gff3
#list of genomes to analyse: genome_files.txt
/data/martin/genomics/analyses/Danaus_genome/Dplex_V4_HiC/braker_annotation_RDK/danaus_plexv4_braker_a006.sequences.tidy.gff3
/data/martin/genomics/analyses/Danaus_genome/Dplex_V4_HiC/Dplex.introns.tidy.gff3
/data/martin/genomics/analyses/Danaus_genome/Dplex_mex/braker_annotation_RDK/danaus_plex_mex_braker_a002.sequences.tidy.gff3
/data/martin/genomics/analyses/Danaus_genome/Dplex_mex/dplex_mex.introns.tidy.gff3
/data/martin/genomics/analyses/Danaus_genome/Dchry2/braker_annotation_RDK_Dchry2_2/Dchry2.2.fa.masked_annotation_a001.sequences.tidy.gff3
#genome_abbs.txt
dplex_v4_B
dplex_v4_O
dplex_mex_B
dplex_mex_O
Dchry2_2
#want to 1)extract exons 2)extract introns 3)calculate line counts of exons 4)count lines of introns 5) get mean exon length 6) get mean intron length
#extract_intron_exon_info.sh
#!/bin/bash
out_dir=/data/martin/genomics/analyses/Danaus_genome/introns_exons/output
echo ${out_dir}
for i in {1..5}
do
# define jobindex
genome_numb=${i}
genome_path=$(cat /data/martin/genomics/analyses/Danaus_genome/introns_exons/genome_files.txt | sed -n ${i}p)
genome_name=$(cat /data/martin/genomics/analyses/Danaus_genome/introns_exons/genome_abbs.txt | sed -n ${i}p)
echo ${genome_path}
echo ${genome_name}
#extract introns/exons
awk '{ if ($3 == "intron") { print } }' ${genome_path} > ${out_dir}/${genome_name}.introns.gff3
awk '{ if ($3 == "exon") { print } }' ${genome_path} > ${out_dir}/${genome_name}.exons.gff3
#get lengths of them
for j in ${out_dir}/${genome_name}.introns.gff3; do
awk '{ print $5-$4; }' "$j" > ${out_dir}/${genome_name}.introns_length_output.txt ;
done
for k in ${out_dir}/${genome_name}.exons.gff3; do
awk '{ print $5-$4; }' "$k" > ${out_dir}/${genome_name}.exons_length_output.txt ;
done
#get means of each
echo ${genome_name} > ${out_dir}/${genome_name}.intron_means.txt
awk '{ total += $1 } END { print total/NR }' ${out_dir}/${genome_name}.introns_length_output.txt >> ${out_dir}/${genome_name}.intron_means.txt
echo ${genome_name} > ${out_dir}/${genome_name}.exon_means.txt
awk '{ total += $1 } END { print total/NR }' ${out_dir}/${genome_name}.exons_length_output.txt >> ${out_dir}/${genome_name}.exon_means.txt
done
cat ${out_dir}/*.intron_means.txt > ${out_dir}/all_introns.txt
cat ${out_dir}/*.exon_means.txt > ${out_dir}/all_exons.txt
####
#out
####
rsync qmaster:/data/martin/genomics/analyses/Danaus_genome/introns_exons/output/*_length_output.txt Dropbox/RishiMAC/Danaus/Genome_annotation/introns/auto_intron_exon/
#cat all_exons.txt
Dchry2_2
217.307
dplex_mex_B
223.69
dplex_mex_O
271.062
dplex_v4_B
208.024
dplex_v4_O
206.286
#cat all_introns.txt
Dchry2_2
1032.27
dplex_mex_B
725.857
dplex_mex_O
559.891
dplex_v4_B
780.619
dplex_v4_O
804.666
#now just for our re-annotated assemblies to allow a direct comparison and only using '.t1' genes
#extract_intron_exon_info_t1.sh
#!/bin/bash
out_dir=/data/martin/genomics/analyses/Danaus_genome/introns_exons/t1_only/output
echo ${out_dir}
for i in {1..3}
do
# define jobindex
genome_numb=${i}
genome_path=$(cat /data/martin/genomics/analyses/Danaus_genome/introns_exons/t1_only/genome_files.txt | sed -n ${i}p)
genome_name=$(cat /data/martin/genomics/analyses/Danaus_genome/introns_exons/t1_only/genome_abbs.txt | sed -n ${i}p)
echo ${genome_path}
echo ${genome_name}
#extract introns/exons
awk '{ if ($3 == "intron") { print } }' ${genome_path} | grep ".t1" > ${out_dir}/${genome_name}.introns.gff3
awk '{ if ($3 == "exon") { print } }' ${genome_path} | grep ".t1" > ${out_dir}/${genome_name}.exons.gff3
#get lengths of them
for j in ${out_dir}/${genome_name}.introns.gff3; do
awk '{ print $5-$4; }' "$j" > ${out_dir}/${genome_name}.introns_length_output.txt ;
done
for k in ${out_dir}/${genome_name}.exons.gff3; do
awk '{ print $5-$4; }' "$k" > ${out_dir}/${genome_name}.exons_length_output.txt ;
done
#get means of each
echo ${genome_name} > ${out_dir}/${genome_name}.intron_means.txt
awk '{ total += $1 } END { print total/NR }' ${out_dir}/${genome_name}.introns_length_output.txt >> ${out_dir}/${genome_name}.intron_means.txt
echo ${genome_name} > ${out_dir}/${genome_name}.exon_means.txt
awk '{ total += $1 } END { print total/NR }' ${out_dir}/${genome_name}.exons_length_output.txt >> ${out_dir}/${genome_name}.exon_means.txt
done
cat ${out_dir}/*.intron_means.txt > ${out_dir}/all_introns.txt
cat ${out_dir}/*.exon_means.txt > ${out_dir}/all_exons.txt
####
#out
####
rsync qmaster:/data/martin/genomics/analyses/Danaus_genome/introns_exons/t1_only/output/*_length_output.txt ./t1_only
#cat all_exons.txt
Dchry2_2
226.447
dplex_mex_B
238.272
dplex_v4_B
216.605
#cat all_introns.txt
Dchry2_2
974.9
dplex_mex_B
665.131
dplex_v4_B
737.858
#######
#functional annotation with pannzer2:
#http://ekhidna2.biocenter.helsinki.fi/barcosel/tmp//RD8AUok6TJT/index.html
#Title: dchrys2.2
#Proteins: 19639
#Database:
#uniprot.Jun2021 consisting of 75694734762 letters and 219740214 sequences
#URL: http://ekhidna2.biocenter.helsinki.fi/barcosel/tmp//RD8AUok6TJT
#Checksum: 265bc5d74e09153deda91b41c5e2c45a7c2529193185fbf98ea5c60e
#how many unique genes are annotated
cut -f1 GO.out | grep ".t1" | sort | uniq | wc -l
#9567
#out of 16260 genes
###########################################################
# 5. BUSCOs:
###########################################################
#busco check on whole assembly and then on proteins:
#/data/martin/genomics/analyses/Danaus_genome/Dchry2/braker_annotation_RDK_Dchry2_2/buscos_G3_paper
conda create -n busco
sconda busco
conda install -c bioconda busco
conda install -c conda-forge parallel
#want buscos for all assemblies and all annotations:
#busco_input_assemblies.txt
/data/martin/genomics/analyses/Danaus_genome/Dchry2/Dchry2.2.fa.masked
/data/martin/genomics/analyses/Danaus_genome/Dplex_V4_HiC/Dapl_Zhan_v3_HiC.RN.uppercase.fasta.renamed.masked
/data/martin/genomics/analyses/Danaus_genome/Dplex_mex/dplex_mex.uppercase.fa.masked
##and all protein sets:
#busco_input_proteins.txt
/data/martin/genomics/analyses/Danaus_genome/Dchry2/braker_annotation_RDK_Dchry2_2/Dchry2.2.fa.masked_annotation_a001.sequences.gt.aa.fasta
/data/martin/genomics/analyses/Danaus_genome/Dplex_mex/dplex_mex.gt.aa.fasta
/data/martin/genomics/analyses/Danaus_genome/Dplex_mex/braker_annotation_RDK/danaus_plex_mex_braker_a002.sequences.gt.aa.fasta
/data/martin/genomics/analyses/Danaus_genome/Dplex_V4_HiC/Dplex.tidy.gff3.gt.aa.fasta
/data/martin/genomics/analyses/Danaus_genome/Dplex_V4_HiC/braker_annotation_RDK/danaus_plexv4_braker_a006.sequences.gt.aa.fasta
echo 'busco -m genome -l insecta -c 1 -i /data/martin/genomics/analyses/Danaus_genome/Dchry2/Dchry2.2.fa.masked -o Dchry2.2.fa.masked.busco_v5.insecta' > busco_g3.command1.txt
parallel -j 1 'qsub -cwd -N busco -V -pe smp 20 -b yes {}' :::: busco_g3.command1.txt
echo 'busco -m genome -l insecta -c 1 -i /data/martin/genomics/analyses/Danaus_genome/Dplex_V4_HiC/Dapl_Zhan_v3_HiC.RN.uppercase.fasta.renamed.masked -o Dapl_Zhan_v3_HiC.RN.uppercase.fasta.renamed.masked.busco_v5.insecta' > busco_g3.command2.txt
parallel -j 1 'qsub -cwd -N busco -V -pe smp 20 -b yes {}' :::: busco_g3.command2.txt
echo 'busco -m genome -l insecta -c 1 -i /data/martin/genomics/analyses/Danaus_genome/Dplex_mex/dplex_mex.uppercase.fa.masked -o test_dplex_mex.uppercase.fa.masked.busco_v5.insecta' > busco_g3.command3.txt
parallel -j 1 'qsub -cwd -N busco -V -pe smp 20 -l h=c2 -b yes {}' :::: busco_g3.command3.txt
echo 'busco -m protein -l insecta -c 1 -i /data/martin/genomics/analyses/Danaus_genome/Dchry2/braker_annotation_RDK_Dchry2_2/Dchry2.2.fa.masked_annotation_a001.sequences.gt.aa.fasta -o Dchry2.2.fa.masked_annotation_a001.sequences.gt.aa.fasta.busco_v5.insecta' > busco_g3.command4.txt
parallel -j 1 'qsub -cwd -N busco -V -pe smp 20 -l h=c2 -b yes {}' :::: busco_g3.command4.txt
echo 'busco -m protein -l insecta -c 1 -i /data/martin/genomics/analyses/Danaus_genome/Dplex_mex/dplex_mex.gt.aa.fasta -o dplex_mex.gt.aa.fasta.busco_v5.insecta' > busco_g3.command5.txt
parallel -j 1 'qsub -cwd -N busco -V -pe smp 20 -l h=c5 -b yes {}' :::: busco_g3.command5.txt
echo 'busco -m protein -l insecta -c 1 -i /data/martin/genomics/analyses/Danaus_genome/Dplex_mex/braker_annotation_RDK/danaus_plex_mex_braker_a002.sequences.gt.aa.fasta -o danaus_plex_mex_braker_a002.sequences.gt.aa.fasta.busco_v5.insecta' > busco_g3.command6.txt
parallel -j 1 'qsub -cwd -N busco -V -pe smp 20 -l h=c5 -b yes {}' :::: busco_g3.command6.txt
echo 'busco -m protein -l insecta -c 1 -i /data/martin/genomics/analyses/Danaus_genome/Dplex_V4_HiC/Dplex.tidy.gff3.gt.aa.fasta -o Dplex.tidy.gff3.gt.aa.fasta.busco_v5.insecta' > busco_g3.command7.txt
parallel -j 1 'qsub -cwd -N busco -V -pe smp 20 -l h=c5 -b yes {}' :::: busco_g3.command7.txt
echo 'busco -m protein -l insecta -c 1 -i /data/martin/genomics/analyses/Danaus_genome/Dplex_V4_HiC/braker_annotation_RDK/danaus_plexv4_braker_a006.sequences.gt.aa.fasta -o danaus_plexv4_braker_a006.sequences.gt.aa.fasta.busco_v5.insecta' > busco_g3.command8.txt
parallel -j 1 'qsub -cwd -N busco -V -pe smp 20 -l h=c2 -b yes {}' :::: busco_g3.command8.txt
for file in *.o*
do
echo ${file}
echo 'busco - output'
grep -A11 "Results from dataset insecta_odb10" ${file}
done
#now test using only the longest transcripts
#/data/martin/genomics/analyses/Danaus_genome/Dchry2/braker_annotation_RDK_Dchry2_2/buscos_G3_paper/longest_transcripts
grep -A1 ".t1" /data/martin/genomics/analyses/Danaus_genome/Dchry2/braker_annotation_RDK_Dchry2_2/Dchry2.2.fa.masked_annotation_a001.sequences.gt.aa.fasta > ./LONGESTDchry2.2.fa.masked_annotation_a001.sequences.gt.aa.fasta
grep -A1 ".t1" /data/martin/genomics/analyses/Danaus_genome/Dplex_mex/braker_annotation_RDK/danaus_plex_mex_braker_a002.sequences.gt.aa.fasta > ./LONGESTdanaus_plex_mex_braker_a002.sequences.gt.aa.fasta
grep -A1 ".t1" /data/martin/genomics/analyses/Danaus_genome/Dplex_V4_HiC/braker_annotation_RDK/danaus_plexv4_braker_a006.sequences.gt.aa.fasta > ./LONGESTdanaus_plexv4_braker_a006.sequences.gt.aa.fasta
echo 'busco -m protein -l insecta -c 1 -i ./LONGESTDchry2.2.fa.masked_annotation_a001.sequences.gt.aa.fasta -o longest_Dchr2.2.busco_v5.insecta' > busco_g3.command9.txt
parallel -j 1 'qsub -cwd -N busco -V -pe smp 20 -l h=c2 -b yes {}' :::: busco_g3.command9.txt
echo 'busco -m protein -l insecta -c 1 -i ./LONGESTdanaus_plex_mex_braker_a002.sequences.gt.aa.fasta -o longest_Dplex_mex.busco_v5.insecta' > busco_g3.command10.txt
parallel -j 1 'qsub -cwd -N busco -V -pe smp 20 -l h=c2 -b yes {}' :::: busco_g3.command10.txt
echo 'busco -m protein -l insecta -c 1 -i ./LONGESTdanaus_plexv4_braker_a006.sequences.gt.aa.fasta -o longest_Dplex_v4.busco_v5.insecta' > busco_g3.command11.txt
parallel -j 1 'qsub -cwd -N busco -V -pe smp 20 -l h=c2 -b yes {}' :::: busco_g3.command11.txt
#checking qv score with merqury
sconda /ceph/users/amackintosh/.conda/envs/assembly
seqtk mergepe /data/martin/genomics/analyses/Danaus_genome/Dchry2/ENA_upload/raw_reads/Illumina/MB18111_Illumina_Dchry2.2_read1.fastq.gz /data/martin/genomics/analyses/Danaus_genome/Dchry2/ENA_upload/raw_reads/Illumina/MB18111_Illumina_Dchry2.2_read2.fastq.gz | gzip > merged_interleaved_MB18111.fq.gz
sconda ceph/users/amackintosh/.conda/envs/merqury/
meryl count k=21 memory=100 threads=16 output reads.meryl ./merged_interleaved_MB18111.fq.gz
$MERQURY/merqury.sh reads.meryl /data/martin/genomics/analyses/Danaus_genome/Dchry2/Dchry2.2.fa Dchry2.2_QV > Dchry2.2_QV.log