Project1-de-novo-assembly_srcipt/Annotation_Diamond at main · Justinma-web/Project1-de-novo-assembly_srcipt · GitHub

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
#!/bin/bash

###1. Annotation using SwissProt
## Download Swissprot database
wget https://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz
unzip uniprot_sprot.fasta.gz

## Build DIAMOND database
diamond makedb --in uniprot_sprot.fasta --db swissprot.dmnd

## Align DNA query sequences against a protein reference database
diamond blastx \
    -d swissprot.dmnd \
    -q cluster_transcripts_fpkm0.5_s3.fa \
    -o sand_smelt_vs_swissprot.xml \
    --threads 20 \
    --max-target-seqs 5 \
    --ultra-sensitive \
    --evalue 1e-5 \
    --outfmt 5

###2. Extract unannotated sequence (Use "blastxml2tab" in this git)
## Filter annotated sequences
wget https://github.com/ISUgenomics/common_scripts/archive/refs/heads/master.zip
unzip master.zip
python2 /home/jade/common_scripts-master/blastXML2Tab.py \
    -o sand_smelt_vs_swissprot.m8 \
    -c std \
    sand_smelt_vs_swissprot.xml

cut -f1 sand_smelt_vs_swissprot.m8 | sort -u > annotated_ids_sws.txt
seqkit grep -v -f annotated_ids_sws.txt cluster_transcripts_fpkm0.5_s3.fa > unannotated_1.fa

##3. Annotate the rest but using TrEMBL
wget https://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_trembl.fasta.gz
gunzip uniprot_trembl.fasta.gz

## Build DIAMOND database
diamond makedb --in uniprot_trembl.fasta --db trembl --threads 20

## Align DNA query sequences against a protein reference database
/home/jade/diamond blastx \
    -d trembl.dmnd \
    -q unannotated_1.fa \
    -o sand_smelt_vs_trembl.xml \
    --threads 20 \
    --max-target-seqs 5 \
    --ultra-sensitive \
    --evalue 1e-5 \
    --outfmt 5

###4. Extract unannotated sequence (Use "blastxml2tab" in this git)
## Filter annotated sequences
python2 /home/jade/common_scripts-master/blastXML2Tab.py \
    -o sand_smelt_vs_trembl.m8 \
    -c std \
    sand_smelt_vs_trembl.xml

cut -f1 sand_smelt_vs_trembl.m8 | sort -u > annotated_ids_tre.txt
seqkit grep -v -f annotated_ids_tre.txt unannotated_1.fa > unannotated_2.fa

##5. Annotate the rest but using uniRef90
wget https://ftp.uniprot.org/pub/databases/uniprot/uniref/uniref90/uniref90.fasta.gz
gunzip uniref90.fasta.gz

## Build DIAMOND database
/home/jade/diamond makedb --in uniref90.fasta.gz --db uniref90 --threads 20

## Align DNA query sequences against a protein reference database
/home/jade/diamond blastx \
    -d uniref90.dmnd \
    -q unannotated_2.fa \
    -o sand_smelt_vs_uniref90.xml \
    --threads 20 \
    --max-target-seqs 5 \
    --ultra-sensitive \
    --evalue 1e-5 \
    --outfmt 5

###5. Extract unannotated sequence (Use "blastxml2tab" in this git)
## Filter annotated sequences
python2 /home/jade/common_scripts-master/blastXML2Tab.py \
    -o sand_smelt_vs_uniref90.m8 \
    -c std \
    sand_smelt_vs_uniref90.xml

cut -f1 sand_smelt_vs_uniref90.m8 | sort -u > annotated_ids_ur90.txt
seqkit grep -v -f annotated_ids_ur90.txt unannotated_2.fa > unannotated_3.fa

##5. Annotate the rest using nr
mkdir -p nr_files
cd nr_files

# Download all nr.*.tar.gz files in parallel with wget
seq -w 000 128 | sed 's|^|https://ftp.ncbi.nlm.nih.gov/blast/db/nr.|;s|$|.tar.gz|' > nr_urls.txt
cat nr_urls.txt | parallel -j 20 "wget -c {}"

# Extract all tarballs in parallel
ls nr.*.tar.gz | parallel -j 20 'tar -xzf {}'

blastdbcmd -db nr -entry all -out nr.fasta -dbtype prot

## Build DIAMOND database
/home/jade/diamond makedb --in nr.fasta --db nr --threads 20

## Align DNA query sequences against a protein reference database
cd /home/jade/1_schooling_transcriptome_annotation
/home/jade/diamond blastx \
    -d /home/jade/1_schooling_transcriptome_annotation/nr_files/nr.dmnd \
    -q unannotated_3.fa \
    -o sand_smelt_vs_nr.xml \
    --threads 20 \
    --max-target-seqs 5 \
    --ultra-sensitive \
    --evalue 1e-5 \
    --outfmt 5

## Integrate all annotated sequence
xxxxxx