@@ -0,0 +1,2609 @@
+# IdentifyOffTargetSiteSequences.py%0A# Shengdar Tsai (stsai4@mgh.harvard.edu)%0A%0A# A program to identify Cas9 off-target sites from in vitro GUIDE-seq data%0A%0A__author__ = 'shengdar'%0A%0Aimport argparse%0Aimport HTSeq%0Aimport itertools%0A# import pyfaidx%0A%0Aparser = argparse.ArgumentParser(description='Identify off-target candidates from Illumina short read sequencing data.')%0Aparser.add_argument('--bam', help='Sorted BAM file', required=True)%0Aargs = parser.parse_args()%0A%0A### Tabulate the start positions for the 2nd read in pair across the genome.%0A### Filter mapped reads with a quality score of 50.%0A%0Adef tabulate_start_positions(BamFileName):%0A    sorted_bam_file = HTSeq.BAM_Reader(args.bam)%0A    ga = HTSeq.GenomicArray(%22auto%22, stranded=False)%0A    ga_windows = HTSeq.GenomicArray(%22auto%22, stranded=False)%0A    ga_stranded = HTSeq.GenomicArray(%22auto%22, stranded=True)%0A%0A    for read in itertools.islice( sorted_bam_file, 1000000 ):  # printing first 5 reads%0A    # for read in sorted_bam_file:%0A        if read.pe_which == %22second%22 and read.aligned and read.aQual %3E= 50:%0A            iv = read.iv%0A            chr = iv.chrom%0A            position = iv.start_d%0A            strand = iv.strand%0A            ga%5BHTSeq.GenomicPosition(chr, position, strand)%5D += 1%0A            ga_windows%5BHTSeq.GenomicPosition(chr, position, strand)%5D = 1%0A            ga_stranded%5BHTSeq.GenomicPosition(chr, position, strand)%5D += 1%0A    return ga, ga_windows, ga_stranded%0A%0A### find_windows method%0Adef find_windows(ga, ga_windows, window_size=3):%0A    # Initialize comparison position%0A    last = HTSeq.GenomicInterval(%220%22, 0, 0)%0A    # Iterate through window GenomicArray and consolidate windows that are within 3 bp%0A    for iv, value in ga_windows.steps():%0A        if value:%0A            if iv.chrom != last.chrom or iv.start - last.end %3E window_size:%0A                last = iv%0A            else:%0A                consolidated_interval = HTSeq.GenomicInterval(iv.chrom, last.start, iv.end)%0A                ga_windows%5Bconsolidated_interval%5D = 1%0A                last = consolidated_interval%0A    # Return consolidated GenomicArray%0A%0A    for iv, value in ga_windows.steps():%0A        if value:%0A            count = sum(list(ga%5Biv%5D))%0A            if count %3E= 3:%0A                print iv, sum(list(ga%5Biv%5D))%0A%0A    return ga_windows%0A%0Adef main():%0A    # Tabulate start positions for read 2%0A    # Identify positions where there are more than 1 read%0A%0A    ga, ga_windows, ga_stranded = tabulate_start_positions(%22/Users/shengdar/Local-projects/invitro_GUIDE-Seq/BAM/SQT01_S1.bam%22)%0A    ga_consolidated_windows = find_windows(ga, ga_windows)%0A%0Aif __name__ == %22__main__%22:%0A    main()%0A%0A%0A

Add FindCleavageSites.py