-
Notifications
You must be signed in to change notification settings - Fork 1
/
depersonalizebam
executable file
·62 lines (60 loc) · 1.91 KB
/
depersonalizebam
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
#!/bin/bash
#
# Depersonalize a bamfile (Remove CIGAR, SEQ, QUAL and additional columns information)
#
# See Samtools Specification: samtools.sourceforge.net/SAMv1.pdf
#
# Usage: ./depersonalizebam input.bam ReadLengthINT > output.bam
#
# 2014-01-30 [email protected]
BAM=$1
READLENGTH=${2:-50}
samtools view -h $BAM | \
awk -v CIGAR=$READLENGTH 'BEGIN \
{
FS="\t"
OFS="\t"
}
{
if (/^@/)
{ # Print header lines without modification
print;
}
else
{
# Update PNEXT (sense strand) or POS (minus strand) from CIGAR
POS=$4
CIG=CIGAR"M"
PNEXT=$8
TLEN=$9
if ($6 !~ /CIG/ && TLEN != 0 && PNEXT != POS)
{ #check if TLEN information is given, then it is paired end
if ( and($2, 0x10) != 16)
{ # sense strand
PNEXT=POS + TLEN - CIGAR
}
else
{ # antisense strand (FAKE START introduction)
POS=PNEXT - TLEN - CIGAR
}
}
if ($NF ~ /YT:Z/)
{
# Put out bowtie alignment score
# YT:Z:<S> Value of UU indicates the read was not part of a pair.
# Value of CP indicates the read was part of a pair and the pair
# aligned concordantly. Value of DP indicates the read was part
# of a pair and the pair aligned discordantly. Value of UP
# indicates the read was part of a pair but the pair failed to
# aligned either concordantly or discordantly.
# (http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml#sam-output)
# Output 11 columns: QNAME FLAG RNAME POS MAPQ CIGAR RNEXT PNEXT TLEN SEQ QUAL YT/YZ
print $1, $2, $3, POS, $5, CIG, "=", PNEXT,TLEN, "*", "*", $NF
}
else
{
# Output 10 columns: QNAME FLAG RNAME POS MAPQ CIGAR RNEXT PNEXT TLEN SEQ QUAL
print $1, $2, $3, POS, $5, CIG, "=", PNEXT,TLEN, "*", "*"
}
}
}' | samtools view -bS -