ProkTNseq - marcottelab/HTseq-toolbox GitHub Wiki

Table of Contents Mapping single-end reads with 'bwa mem' or 'bwa sw' Filter best hits Make tnHit and gHit files Calculate tnSite frequency

• https://github.com/marcottelab/HTseq-toolbox/blob/master/sam/sam-to-sam_hit.py

#!/bin/bash
BWA="$HOME/src.HTseq/bwa/0.7.10/bwa"
SAM2HIT="$HOME/git/HTseq-toolbox/sam/sam-to-sam_hit.py"

THREADS=12

DB="$WORK/PSEAE.db/bwadb/PSEAE_PA14_genome"
DBNAME=$(basename $DB)
DBNAME=${DBNAME/.fa/}

for FQ in $(ls ../fastq.tn/*fastq)
do
  SAM=$(basename $FQ)
  SAM=${SAM/.called.fastq/_called}
  SAM=$SAM"."$DBNAME".bwa_sw.sam"

  $BWA bwasw -t $THREADS $DB $FQ > $SAM
  $SAM2HIT $SAM
done

• https://github.com/marcottelab/HTseq-toolbox/blob/master/sam/sam-to-sam_best.py

#!/bin/bash
SAM2BEST="$HOME/git/HTseq-toolbox/sam/sam-to-sam_best.py"

for SAM in $(ls *sam_hit)
do
  $SAM2BEST $SAM
done

• https://github.com/marcottelab/HTseq-toolbox/blob/master/sam/tnseq_sam-to-tnHit+gHit.py

#!/bin/bash
TNHIT="$HOME/git/HTseq-toolbox/sam/tnseq_sam-to-tnHit+gHit.py"

for SAM in $(ls *sam_hit_best)
do
  $TNHIT $SAM
done

• Output gHit file format (reads almost perfectly matched to the genome)

HWI-D00289:132:C4Y44ACXX:1:1104:4287:25326	-	5739694	5739776	CTGTCGCCCTGCGACCTTCTGCGCATTGGCGGGCGCCGCCACCAGCGCGTAGGGTAACGCGCTCTCTCTCCTCGGCCACAAGTAGCGTCCTGGACGAGCAG
HWI-D00289:132:C4Y44ACXX:1:1105:14692:95713	+	5688961	5689043	TTAAACGTCCAGGACGCTACTTGTGCGAAGGGTTAGGCAGGAGAAAGCGTCAGCCGCGCATCGGGCGGCGACCGCTCTCAGGGTTCTCTTCGATCTGCGTC

• Output tnHit file format (reads partially matched to the genome; uppercase bases are the part of a read matched to the genome, and the lowercase bases are not matched to the genome)

HWI-D00289:132:C4Y44ACXX:1:1101:10000:91751	-	2777872	2777937	GGCGATCCTGGTGTCCCTGGTGGTATCCCTGACCCTCACGCCGATGCTCTGCGCGCGTCTGCTGCtgactcttatacacaagtagcgtcctggacgctgaa
HWI-D00289:132:C4Y44ACXX:1:1101:10001:25138	-	5777862	5777916	ggggggggccCCCGGGTCTTCGCCGAAAGCTCGCCGGCCAGGGCCTGCAGACGCTCCTCGCGACctgactcttatacacaagtagcgtcctggacgttaga
HWI-D00289:132:C4Y44ACXX:1:1101:10002:68489	+	4170170	4170234	ctctgcgtccaggacgctacttgtgtataagagtcagCCTCCAGGGGCGCGCGCGATACCCCCACCAGGCTGACGTCGTGGGGCGCCCCGTGCCGGCTGGC

• https://github.com/marcottelab/HTseq-toolbox/blob/master/sam/tnHit-to-tnSite.py

 #!/bin/bash
for HIT in $(ls *tnHit)
do
  echo $HIT
  $HOME/git/HTseq-toolbox/sam/tnHit-to-tnSite.py $HIT taagagtcag
done

• Output tnSite file format (genomic position - forward read count - reverse read count)

185	150	0
217	0	38
218	21	0
228	0	459
230	2	0
257	0	5

• Output tnSite.log file format

# TN_f: taagagtcag
# TN_r: ctgactctta
#max_mismatch=1
HWI-D00289:132:C4Y44ACXX:1:1101:10001:40772	adjust_1	+	ataagagtca.GCTCTTGTGC -> taagagtcag.CTCTTGTGC
HWI-D00289:132:C4Y44ACXX:1:1101:10004:77184	adjust_1	+	ataagagtca.GGAGCAGGAG -> taagagtcag.GAGCAGGAG
....
HWI-D00289:132:C4Y44ACXX:8:2309:9999:37562	adjust_1	-	GCTCTACGCC.tgactcttat -> GCTCTACGC.ctgactctta
# Total hits: 22235051
# Perfect edges: 12412119 (55.82 pct)
# Substituted edges: 210726 (0.95 pct)
# Adjusted edges: 6239315 (28.06 pct)
# Errors: 3372891 (15.17 pct)
# Unique hits: 19098516