-
Notifications
You must be signed in to change notification settings - Fork 9
/
hugeseq.sh
executable file
·94 lines (82 loc) · 3.87 KB
/
hugeseq.sh
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
#!/bin/sh -l
ref_dir=/srv/gs1/projects/cherry/sacCer3
name=$1
seq_dir=$2
out_dir=$3
if [ $# -ne 3 ] && [ $# -ne 4 ]
then
echo "Usage: $SCRIPT strain_name seq_dir output_dir [32 or 64]"
exit 1
elif [ $# == 4 ]
then
qual_encode=$4
else
qual_encode=32
fi
if [ $qual_encode -ne 32 ] && [ $qual_encode -ne 64 ]
then
echo "Usage: $SCRIPT strain_name seq_dir output_dir [32 or 64]"
exit 1
fi
prep.sh $seq_dir/"$name"_1.fastq $out_dir/"$name"_1.fastq
prep.sh $seq_dir/"$name"_2.fastq $out_dir/"$name"_2.fastq
if [ $qual_encode == 64 ]
then
bwa aln -q 15 -l 35 -k 2 -n 0.04 -o 2 -e 6 -t 1 -I $ref_dir/sacCer3.fa "$seq_dir/$name"_1.fastq > "$out_dir/$name"_1.fastq.sai
bwa aln -q 15 -l 35 -k 2 -n 0.04 -o 2 -e 6 -t 1 -I $ref_dir/sacCer3.fa "$seq_dir/$name"_2.fastq > "$out_dir/$name"_2.fastq.sai
else
bwa aln -q 15 -l 35 -k 2 -n 0.04 -o 2 -e 6 -t 1 $ref_dir/sacCer3.fa "$seq_dir/$name"_1.fastq > "$out_dir/$name"_1.fastq.sai
bwa aln -q 15 -l 35 -k 2 -n 0.04 -o 2 -e 6 -t 1 $ref_dir/sacCer3.fa "$seq_dir/$name"_2.fastq > "$out_dir/$name"_2.fastq.sai
fi
aln_sam.sh $out_dir/$name.bam "$out_dir/$name"_1.fastq "$out_dir/$name"_2.fastq "@RG\tID:ILLUMINA_"$name"\tLB:3\tSM:0167\tPL:illumina"
sam_sort.sh $out_dir/$name.bam $out_dir/$name.sorted.bam 12
sam_index.sh $out_dir/$name.sorted.bam
sam_rm.sh $out_dir/$name.bam
while read line
do
chr=`echo $line | awk '{print $1}'`
bin_sam.sh $chr $out_dir/$chr.bam $out_dir/$name.sorted.bam
sam_index.sh $out_dir/$chr.bam
clean_nodup.sh $out_dir/$chr.bam $out_dir/$chr.nodup.bam
sam_index.sh $out_dir/$chr.nodup.bam
rm -rf $out_dir/$chr.bam $out_dir/$chr.bam.bai
clean_realn.sh $out_dir/$chr.nodup.bam $out_dir/$chr.realn.bam
sam_index.sh $out_dir/$chr.realn.bam
rm -rf $out_dir/$chr.nodup.bam $out_dir/$chr.nodup.bam.bai
clean_recal.sh $out_dir/$chr.realn.bam $out_dir/$chr.recal.temp.bam
java -Xms5g -Xmx5g -jar $PICARD/AddOrReplaceReadGroups.jar INPUT=$out_dir/$chr.recal.temp.bam OUTPUT=$out_dir/$chr.recal.bam RGID=ILLUMINA_$name RGPL=illumina RGLB=1 RGPU=1111 RGSM=$name VALIDATION_STRINGENCY=LENIENT
sam_index.sh $out_dir/$chr.recal.bam
rm -rf $out_dir/$chr.realn.bam $out_dir/$chr.realn.bam.bai $out_dir/$chr.realn.intervals $out_dir/$chr.recal.temp.bam
done < $ref_dir/sacCer3.fa.fai
snp_vcfs=''
indel_vcfs=''
pileup_vcfs=''
cat $ref_dir/sacCer3.fa.fai | ( while read line
do
chr=`echo $line | awk '{print $1}'`
var_snp.sh $out_dir/$chr.snp.gatk.vcf $out_dir/$chr.recal.bam
var_filter.sh $out_dir/$chr.snp.gatk.vcf
var_indel.sh $out_dir/$chr.indel.gatk.vcf $out_dir/$chr.recal.bam
var_filter.sh $out_dir/$chr.indel.gatk.vcf
var_pileup.sh $out_dir/$chr.pileup.vcf $out_dir/$chr.recal.bam
snp_vcfs="$snp_vcfs $out_dir/$chr.snp.gatk.vcf"
indel_vcfs="$indel_vcfs $out_dir/$chr.indel.gatk.vcf"
pileup_vcfs="$pileup_vcfs $out_dir/$chr.pileup.gatk.vcf"
done
concat_vcf.sh $out_dir/$name.snp.gatk.raw.vcf $snp_vcfs
concat_vcf.sh $out_dir/$name.indel.gatk.raw.vcf $indel_vcfs
concat_vcf.sh $out_dir/$name.pileup.raw.vcf $pileup_vcfs
#merge_vcf.sh $out_dir/$name.snpindel.vcf $out_dir/$name.snp.gatk.raw.vcf $out_dir/$name.indel.gatk.raw.vcf $out_dir/$name.pileup.raw.vcf
annotate.py $out_dir/$name.snp.tsv $out_dir/$name.snp.gatk.raw.vcf
annotate.py $out_dir/$name.indel.tsv $out_dir/$name.indel.gatk.raw.vcf
annotate.py $out_dir/$name.pileup.tsv $out_dir/$name.pileup.raw.vcf
gffs=''
chr=`echo $line | awk '{print $1}'`
var_sv_rpm.sh $out_dir/$name.rpm.gff $out_dir/$name.sorted.bam
var_sv_sra.sh $out_dir/$name.sra.gff $out_dir/$name.sorted.bam
var_sv_rda.sh $out_dir/$name.rda.gff $out_dir/$name.sorted.bam
var_sv_jct.sh $out_dir/$name.jct.gff $out_dir/$name.sorted.bam
var_sv_ctx.sh $out_dir/$name.ctx.gff $out_dir/$name.sorted.bam
gffs="$gffs $out_dir/$name.rpm.gff $out_dir/$name.sra.gff $out_dir/$name.rda.gff $out_dir/$name.jct.gff $out_dir/$name.ctx.gff"
merge_gff.sh $out_dir/$name.svcnv.gff $gffs
annotate.py $out_dir/$name.svcnv.tsv $out_dir/$name.svcnv.raw.gff )