pyBigWig
0.3.23
ส่วนขยาย Python ที่เขียนด้วย C เพื่อเข้าถึงไฟล์ Bigbed อย่างรวดเร็วและการเข้าถึงและการสร้างไฟล์ BigWig ส่วนขยายนี้ใช้ libbigwig สำหรับการเข้าถึงไฟล์ในพื้นที่และระยะไกล
คุณสามารถติดตั้งส่วนขยายนี้ได้โดยตรงจาก GitHub ด้วย:
pip install pyBigWig
หรือกับ conda
conda install pybigwig -c conda-forge -c bioconda
ต้องติดตั้งข้อกำหนดที่ไม่ใช่ Python ตาม:
curl-config config)จำเป็นต้องมีส่วนหัวและห้องสมุดสำหรับสิ่งเหล่านี้
การใช้งานขั้นพื้นฐานมีดังนี้:
>>> import pyBigWig
สิ่งนี้จะทำงานได้หากไดเรกทอรีการทำงานของคุณคือไดเรกทอรีซอร์สโค้ด Pybigwig
>>> bw = pyBigWig.open("test/test.bw")
โปรดทราบว่าหากไฟล์ไม่มีอยู่คุณจะเห็นข้อความแสดงข้อผิดพลาดและจะ None การส่งคืน เป็นค่าเริ่มต้นไฟล์ทั้งหมดจะเปิดสำหรับการอ่านและไม่เขียน คุณสามารถเปลี่ยนแปลงสิ่งนี้ได้โดยผ่านโหมดที่มี w :
>>> bw = pyBigWig.open("test/output.bw", "w")
โปรดทราบว่าไฟล์ที่เปิดสำหรับการเขียนไม่สามารถสอบถามได้สำหรับช่วงเวลาหรือสถิติของมันสามารถเขียนได้ เท่านั้น หากคุณเปิดไฟล์สำหรับการเขียนคุณจะต้องเพิ่มส่วนหัว (ดูส่วนด้านล่างนี้)
นอกจากนี้ยังรองรับการเข้าถึงการอ่าน Bigbed ในพื้นที่และระยะไกล:
>>> bb = pyBigWig.open("https://www.encodeproject.org/files/ENCFF001JBR/@@download/ENCFF001JBR.bigBed")
ในขณะที่คุณสามารถระบุโหมดสำหรับไฟล์ bigbed ได้ แต่ก็ถูกละเว้น วัตถุที่ส่งคืนโดย pyBigWig.open() จะเหมือนกันไม่ว่าคุณจะเปิดไฟล์ bigwig หรือ bigbed
เนื่องจากทั้งสองไฟล์ Bigwig และ Bigbed สามารถเปิดได้จึงจำเป็นต้องพิจารณาว่าวัตถุ bigWigFile ที่กำหนดไปยังไฟล์ bigwig หรือ bigbed หรือไม่ ด้วยเหตุนี้เราจึงสามารถใช้ฟังก์ชัน isBigWig() และ isBigBed() :
>>> bw = pyBigWig.open("test/test.bw")
>>> bw.isBigWig()
True
>>> bw.isBigBed()
False
วัตถุ bigWigFile มีพจนานุกรมที่ถือความยาวโครโมโซมซึ่งสามารถเข้าถึงได้ด้วย accessor chroms()
>>> bw.chroms()
dict_proxy({'1': 195471971L, '10': 130694993L})
นอกจากนี้คุณยังสามารถสอบถามโครโมโซมเฉพาะโดยตรงได้
>>> bw.chroms("1")
195471971L
ความยาวจะถูกเก็บไว้เป็นจำนวนเต็ม "ยาว" ซึ่งเป็นเหตุผลว่าทำไมจึงมีคำต่อท้าย L หากคุณระบุโครโมโซมที่ไม่มีอยู่จริง
>>> bw.chroms("c")
>>>
บางครั้งมันมีประโยชน์ในการพิมพ์ส่วนหัวของ Bigwig สิ่งนี้ถูกนำเสนอที่นี่เป็นพจนานุกรม Python ที่มี: รุ่น (โดยทั่วไปคือ 4 ), จำนวนของระดับการซูม ( nLevels ), จำนวนฐานที่อธิบาย ( nBasesCovered ), ค่าต่ำสุด ( minVal ), ค่าสูงสุด ( maxVal ), The ผลรวมของค่าทั้งหมด ( sumData ) และผลรวมของค่ากำลังสองทั้งหมด ( sumSquared ) สองสิ่งสุดท้ายนี้จำเป็นสำหรับการพิจารณาค่าเฉลี่ยและค่าเบี่ยงเบนมาตรฐาน
>>> bw.header()
{'maxVal': 2L, 'sumData': 272L, 'minVal': 0L, 'version': 4L, 'sumSquared': 500L, 'nLevels': 1L, 'nBasesCovered': 154L}
โปรดทราบว่าสิ่งนี้เป็นไปได้สำหรับไฟล์ bigbed และคีย์พจนานุกรมเดียวกันจะปรากฏขึ้น รายการต่าง ๆ เช่น maxVal , sumData , minVal และ sumSquared นั้นส่วนใหญ่ไม่มีความหมาย
ไฟล์ BigWig ใช้เพื่อจัดเก็บค่าที่เกี่ยวข้องกับตำแหน่งและช่วงของมัน โดยทั่วไปเราต้องการเข้าถึงค่าเฉลี่ยในช่วงอย่างรวดเร็วซึ่งง่ายมาก:
>>> bw.stats("1", 0, 3)
[0.2000000054637591]
สมมติว่าแทนที่จะเป็นค่าเฉลี่ยเราต้องการค่าสูงสุดแทน:
>>> bw.stats("1", 0, 3, type="max")
[0.30000001192092896]
ตัวเลือกอื่น ๆ คือ "ขั้นต่ำ" (ค่าต่ำสุด), "ความครอบคลุม" (ส่วนของฐานที่ครอบคลุม) และ "std" (ค่าเบี่ยงเบนมาตรฐานของค่า)
บ่อยครั้งที่เราต้องการคำนวณค่าของถังขยะระยะห่างเท่า ๆ กันในช่วงเวลาที่กำหนดซึ่งเป็นเรื่องง่ายเช่นกัน:
>>> bw.stats("1",99, 200, type="max", nBins=2)
[1.399999976158142, 1.5]
nBins เริ่มต้นเป็น 1 เช่นเดียวกับ type เริ่มต้นที่จะ mean
หากไม่ได้รับตำแหน่งเริ่มต้นและสิ้นสุดการใช้โครโมโซมทั้งหมด:
>>> bw.stats("1")
[1.3351851569281683]
หมายเหตุถึงผู้อ่าน Lay: ส่วนนี้ค่อนข้างทางเทคนิคและรวมอยู่ในความสมบูรณ์เท่านั้น สรุปคือหากความต้องการของคุณต้องการค่าเฉลี่ย/สูงสุด/ฯลฯ ค่าสรุปสำหรับช่วงเวลาหรือช่วงเวลาและการแลกเปลี่ยนความเร็วเล็กน้อยเป็นที่ยอมรับได้ว่าคุณควรใช้ตัวเลือก
exact=Trueในฟังก์ชันstats()
โดยค่าเริ่มต้นมีบางแง่มุมที่ไม่ได้ใช้งานง่ายในการคำนวณสถิติเกี่ยวกับช่วงในไฟล์ Bigwig รูปแบบ Bigwig ถูกสร้างขึ้นในบริบทของเบราว์เซอร์จีโนม ที่นั่นการคำนวณสถิติสรุปที่แน่นอนสำหรับช่วงเวลาที่กำหนดมีความสำคัญน้อยกว่าการคำนวณสถิติโดยประมาณอย่างรวดเร็ว (หลังจากทั้งหมดเบราว์เซอร์จะต้องสามารถแสดงช่วงเวลาที่ต่อเนื่องกันและสนับสนุนการเลื่อน/ซูม) ได้อย่างรวดเร็ว ด้วยเหตุนี้ไฟล์ BigWig จึงมีความสัมพันธ์ระหว่างค่าช่วงเวลา แต่ยังรวมถึง sum of values / sum of squared values / minimum value / maximum value / number of bases covered สำหรับถังขยะขนาดเท่ากันในขนาดต่างๆ ขนาดที่แตกต่างกันเหล่านี้เรียกว่า "ระดับการซูม" ระดับการซูมที่เล็กที่สุดมีถังขยะที่มีขนาดเฉลี่ย 16 เท่าในไฟล์และระดับซูมแต่ละระดับที่ตามมามีถังขยะที่ใหญ่กว่าก่อนหน้า 4 เท่า วิธีการนี้ใช้ในเครื่องมือของเคนท์ดังนั้นจึงมีแนวโน้มที่จะใช้ในเกือบทุกไฟล์ BigWig ที่มีอยู่ในปัจจุบัน
เมื่อไฟล์ bigwig ถูกสอบถามสำหรับสถิติสรุปขนาดของช่วงเวลาจะถูกใช้เพื่อตรวจสอบว่าจะใช้ระดับการซูมและถ้าเป็นเช่นนั้น ระดับการซูมที่ดีที่สุดคือสิ่งที่มีถังขยะที่ใหญ่ที่สุดไม่เกินครึ่งหนึ่งของความกว้างของช่วงเวลาที่ต้องการ หากไม่มีระดับการซูมดังกล่าวจะใช้ช่วงเวลาเดิมแทนสำหรับการคำนวณ
เพื่อความสอดคล้องกับเครื่องมืออื่น ๆ Pybigwig ใช้วิธีการเดียวกันนี้ อย่างไรก็ตามเนื่องจากนี่คือ (a) ไม่ได้ใช้งานง่ายและ (b) ที่ไม่พึงประสงค์ในบางแอปพลิเคชัน Pybigwig ช่วยให้การคำนวณสถิติสรุปที่แน่นอนโดยไม่คำนึงถึงขนาดช่วงเวลา (เช่นจะอนุญาตให้ละเว้นระดับการซูม) สิ่งนี้ถูกเสนอเดิมที่นี่และตัวอย่างด้านล่าง:
>>> import pyBigWig
>>> from numpy import mean
>>> bw = pyBigWig.open("http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeMapability/wgEncodeCrgMapabilityAlign75mer.bigWig")
>>> bw.stats('chr1', 89294, 91629)
[0.20120902053804418]
>>> mean(bw.values('chr1', 89294, 91629))
0.22213841940688142
>>> bw.stats('chr1', 89294, 91629, exact=True)
[0.22213841940688142]
ในขณะที่วิธี stats() สามารถ ใช้เพื่อดึงค่าดั้งเดิมสำหรับแต่ละฐาน (เช่นโดยการตั้ง nBins เป็นจำนวนฐาน) มันจะดีกว่าที่จะใช้ accessor values() แทน
>>> bw.values("1", 0, 3)
[0.10000000149011612, 0.20000000298023224, 0.30000001192092896]
รายการที่ผลิตจะมีค่าหนึ่งค่าสำหรับทุกฐานในช่วงที่ระบุ หากฐานเฉพาะไม่มีค่าที่เกี่ยวข้องในไฟล์ BigWig ค่าที่ส่งคืนจะเป็น nan
>>> bw.values("1", 0, 4)
[0.10000000149011612, 0.20000000298023224, 0.30000001192092896, nan]
บางครั้งก็สะดวกในการดึงรายการทั้งหมดที่ทับซ้อนกันบางช่วง สามารถทำได้ด้วยฟังก์ชั่น intervals() :
>>> bw.intervals("1", 0, 3)
((0, 1, 0.10000000149011612), (1, 2, 0.20000000298023224), (2, 3, 0.30000001192092896))
สิ่งที่ส่งคืนคือรายการของ tuples ที่มี: ตำแหน่งเริ่มต้นตำแหน่งสิ้นสุดและค่า ดังนั้นตัวอย่างข้างต้นมีค่า 0.1 , 0.2 และ 0.3 ที่ตำแหน่ง 0 , 1 และ 2 ตามลำดับ
หากไม่ได้รับตำแหน่งเริ่มต้นและสิ้นสุดช่วงเวลาทั้งหมดในโครโมโซมที่ระบุจะถูกส่งคืน:
>>> bw.intervals("1")
((0, 1, 0.10000000149011612), (1, 2, 0.20000000298023224), (2, 3, 0.30000001192092896), (100, 150, 1.399999976158142), (150, 151, 1.5))
ตรงข้ามกับไฟล์ BigWig ไฟล์ BigBed ถือรายการซึ่งเป็นช่วงเวลาที่มีสตริงที่เกี่ยวข้อง คุณสามารถเข้าถึงรายการเหล่านี้ได้โดยใช้ฟังก์ชั่น entries() :
>>> bb = pyBigWig.open("https://www.encodeproject.org/files/ENCFF001JBR/@@download/ENCFF001JBR.bigBed")
>>> bb.entries('chr1', 10000000, 10020000)
[(10009333, 10009640, '61035t130t-t0.026t0.42t404'), (10014007, 10014289, '61047t136t-t0.029t0.42t404'), (10014373, 10024307, '61048t630t-t5.420t0.00t2672399')]
เอาต์พุตเป็นรายการของ tuples รายการ องค์ประกอบ tuple คือตำแหน่ง start และ end ของแต่ละรายการตามด้วย string ที่เกี่ยวข้อง สตริงจะถูกส่งคืนอย่างตรงไปตรงมาตามที่จัดขึ้นในไฟล์ bigbed ดังนั้นการแยกวิเคราะห์มันจะถูกทิ้งไว้ให้คุณ ในการพิจารณาว่าฟิลด์ต่าง ๆ อยู่ในสตริงเหล่านี้อย่างไรให้ปรึกษาสตริง SQL:
>>> bb.SQL()
table RnaElements
"BED6 + 3 scores for RNA Elements data"
(
string chrom; "Reference sequence chromosome or scaffold"
uint chromStart; "Start position in chromosome"
uint chromEnd; "End position in chromosome"
string name; "Name of item"
uint score; "Normalized score from 0-1000"
char[1] strand; "+ or - or . for unknown"
float level; "Expression level such as RPKM or FPKM. Set to -1 for no data."
float signif; "Statistical significance such as IDR. Set to -1 for no data."
uint score2; "Additional measurement/count e.g. number of reads. Set to 0 for no data."
)
โปรดทราบว่ารายการสามรายการแรกในสตริง SQL ไม่ได้เป็นส่วนหนึ่งของสตริง
หากคุณจำเป็นต้องรู้ว่ารายการอยู่ที่ไหนและไม่ใช่ค่าที่เกี่ยวข้องคุณสามารถบันทึกหน่วยความจำได้โดยการระบุ withString=False ใน entries() :
>>> bb.entries('chr1', 10000000, 10020000, withString=False)
[(10009333, 10009640), (10014007, 10014289), (10014373, 10024307)]
หากคุณเปิดไฟล์เพื่อเขียนคุณจะต้องให้ส่วนหัวก่อนที่คุณจะสามารถเพิ่มรายการใด ๆ ส่วนหัวมีโครโมโซมทั้งหมด ตามลำดับ และขนาดของพวกเขา หากจีโนมของคุณมีโครโมโซมสองตัวคือ CHR1 และ CHR2 ความยาว 1 และ 1.5 ล้านฐานจากนั้นต่อไปนี้จะเพิ่มส่วนหัวที่เหมาะสม:
>>> bw.addHeader([("chr1", 1000000), ("chr2", 1500000)])
ส่วนหัวของ Bigwig นั้นมีความอ่อนไหวเป็นกรณี ๆ ดังนั้น chr1 และ Chr1 จึงแตกต่างกัน ในทำนองเดียวกัน 1 และ chr1 ไม่เหมือนกันดังนั้นคุณจึงไม่สามารถผสมชื่อโครโมโซม Ensembl และ UCSC ได้ หลังจากเพิ่มส่วนหัวแล้วคุณสามารถเพิ่มรายการได้
โดยค่าเริ่มต้นจะมีการสร้าง "ระดับการซูม" สูงสุด 10 รายการสำหรับไฟล์ BigWig คุณสามารถเปลี่ยนหมายเลขเริ่มต้นนี้ด้วยอาร์กิวเมนต์ตัวเลือก maxZooms การใช้งานร่วมกันคือการสร้างไฟล์ bigwig ที่เพียงแค่มีช่วงเวลาและไม่มีระดับการซูม:
>>> bw.addHeader([("chr1", 1000000), ("chr2", 1500000)], maxZooms=0)
หากคุณตั้งค่า maxTooms=0 โปรดทราบว่า IGV และเครื่องมืออื่น ๆ อีกมากมายจะไม่ทำงานเนื่องจากพวกเขาคิดว่าอย่างน้อยหนึ่งระดับการซูมจะปรากฏขึ้น คุณควรใช้ค่าเริ่มต้นเว้นแต่คุณจะไม่คาดหวังว่าไฟล์ Bigwig จะใช้งานโดยแพ็คเกจอื่น ๆ
สมมติว่าคุณเปิดไฟล์สำหรับการเขียนและเพิ่มส่วนหัวคุณสามารถเพิ่มรายการได้ โปรดทราบว่า จะต้อง เพิ่มรายการตามลำดับเนื่องจากไฟล์ BigWig จะมีช่วงเวลาที่สั่งซื้อเสมอ มีสามรูปแบบที่ไฟล์ bigwig สามารถใช้ภายในเพื่อจัดเก็บรายการ รูปแบบที่สังเกตได้มากที่สุดนั้นเหมือนกับไฟล์ bedgraph:
chr1 0 100 0.0
chr1 100 120 1.0
chr1 125 126 200.0
รายการเหล่านี้จะถูกเพิ่มดังนี้:
>>> bw.addEntries(["chr1", "chr1", "chr1"], [0, 100, 125], ends=[5, 120, 126], values=[0.0, 1.0, 200.0])
แต่ละรายการมี 12 ไบต์ก่อนการบีบอัด
รูปแบบที่สองใช้ช่วงคงที่ แต่ขนาดขั้นตอนตัวแปรระหว่างรายการ สิ่งเหล่านี้สามารถแสดงได้ในไฟล์กระดิกเป็น:
variableStep chrom=chr1 span=20
500 -2.0
600 150.0
635 25.0
รายการข้างต้นอธิบาย (1-based) ตำแหน่ง 501-520, 601-620 และ 636-655 สิ่งเหล่านี้จะถูกเพิ่มดังนี้:
>>> bw.addEntries("chr1", [500, 600, 635], values=[-2.0, 150.0, 25.0], span=20)
แต่ละรายการของประเภทนี้มี 8 ไบต์ก่อนการบีบอัด
รูปแบบสุดท้ายใช้ขั้นตอนที่คงที่และขยายสำหรับแต่ละรายการซึ่งสอดคล้องกับรูปแบบการงอ StenedStep:
fixedStep chrom=chr1 step=30 span=20
-5.0
-20.0
25.0
รายการข้างต้นอธิบาย (1-based) ฐาน 901-920, 931-950 และ 961-980 และจะเพิ่มดังนี้:
>>> bw.addEntries("chr1", 900, values=[-5.0, -20.0, 25.0], span=20, step=30)
แต่ละรายการของประเภทนี้มี 4 ไบต์
โปรดทราบว่า Pybigwig จะพยายามป้องกันไม่ให้คุณเพิ่มรายการในลำดับที่ไม่ถูกต้อง อย่างไรก็ตามสิ่งนี้ต้องใช้หัวเพิ่มเติม หากไม่สามารถยอมรับได้คุณสามารถระบุ validate=False เมื่อเพิ่มรายการ:
>>> bw.addEntries(["chr1", "chr1", "chr1"], [100, 0, 125], ends=[120, 5, 126], values=[0.0, 1.0, 200.0], validate=False)
เห็นได้ชัดว่าคุณต้องรับผิดชอบในการตรวจสอบให้แน่ใจว่าคุณ ไม่ได้ เพิ่มรายการตามลำดับ ไฟล์ที่ได้จะไม่สามารถใช้งานได้
ไฟล์สามารถปิดได้ด้วย bw.close() อย่างง่ายตามปกติกับไฟล์ประเภทอื่น สำหรับไฟล์ที่เปิดสำหรับการเขียนการปิดไฟล์จะเขียนรายการบัฟเฟอร์ใด ๆ ไปยังดิสก์สร้างและเขียนดัชนีไฟล์และสร้างระดับการซูม ดังนั้นสิ่งนี้อาจใช้เวลาสักครู่
ณ รุ่น 0.3.0, Pybigwig รองรับอินพุตของพิกัดโดยใช้จำนวนเต็มและเวกเตอร์ numpy ในบางฟังก์ชั่น หากติดตั้ง numpy ก่อนที่จะติดตั้ง pybigwig เพื่อตรวจสอบว่ามีการติดตั้ง pybigwig ด้วยการสนับสนุน NUMPY หรือไม่โดยการตรวจสอบ numpy Accessor:
>>> import pyBigWig
>>> pyBigWig.numpy
1
ถ้า pyBigWig.numpy คือ 1 ดังนั้น Pybigwig จะถูกรวบรวมด้วยการสนับสนุน numpy ซึ่งหมายความว่า addEntries() สามารถยอมรับพิกัด numpy:
>>> import pyBigWig
>>> import numpy
>>> bw = pyBigWig.open("/tmp/delete.bw", "w")
>>> bw.addHeader([("1", 1000)], maxZooms=0)
>>> chroms = np.array(["1"] * 10)
>>> starts = np.array([0, 10, 20, 30, 40, 50, 60, 70, 80, 90], dtype=np.int64)
>>> ends = np.array([5, 15, 25, 35, 45, 55, 65, 75, 85, 95], dtype=np.int64)
>>> values0 = np.array(np.random.random_sample(10), dtype=np.float64)
>>> bw.addEntries(chroms, starts, ends=ends, values=values0)
>>> bw.close()
นอกจากนี้ values() สามารถส่งออกเวกเตอร์ numpy โดยตรง:
>>> bw = bw.open("/tmp/delete.bw")
>>> bw.values('1', 0, 10, numpy=True)
[ 0.74336642 0.74336642 0.74336642 0.74336642 0.74336642 nan
nan nan nan nan]
>>> type(bw.values('1', 0, 10, numpy=True))
<type 'numpy.ndarray'>
หากคุณไม่ได้ติดตั้ง Curl Pybigwig จะถูกติดตั้งโดยไม่สามารถเข้าถึงไฟล์ระยะไกลได้ คุณสามารถตรวจสอบได้ว่าคุณจะสามารถเข้าถึงไฟล์ระยะไกลด้วย pyBigWig.remote ได้หรือไม่ หากส่งคืน 1 คุณสามารถเข้าถึงไฟล์ระยะไกล ถ้ามันส่งคืน 0 คุณก็ทำไม่ได้
ตั้งแต่รุ่น 0.3.5 Pybigwig สามารถอ่านและเขียนไฟล์ Bigwig ที่ขาดรายการได้ โปรดทราบว่าโดยทั่วไปแล้วไฟล์ดังกล่าวจะไม่สามารถใช้งานได้กับโปรแกรมอื่น ๆ เนื่องจากไม่มีคำจำกัดความว่าไฟล์ BigWig ที่ไม่มีรายการควรดูอย่างไร สำหรับไฟล์ดังกล่าว Accessor intervals() จะ None ส่งคืนฟังก์ชัน stats() จะส่งคืนรายการของ None ยาวที่ต้องการและ values() จะส่งคืน [] (รายการว่าง) โดยทั่วไปควรอนุญาตให้โปรแกรมที่ใช้ Pybigwig ดำเนินการต่อโดยไม่มีปัญหา
สำหรับผู้ที่ต้องการเลียนแบบการทำงานของ Pybigwig/Libbigwig ในเรื่องนี้โปรดทราบว่ามันดูจำนวนฐานที่ครอบคลุม (ตามที่รายงานในส่วนหัวไฟล์) เพื่อตรวจสอบไฟล์ "ว่าง"
ไฟล์ Wiggle, Bigwig และ Bigbed ใช้พิกัดครึ่งเปิดที่ใช้ 0 ซึ่งใช้โดยส่วนขยายนี้ ดังนั้นในการเข้าถึงค่าสำหรับฐานแรกของ chr1 หนึ่งจะระบุตำแหน่งเริ่มต้นเป็น 0 และตำแหน่งสิ้นสุดเป็น 1 ในทำนองเดียวกันฐาน 100 ถึง 115 จะมีจุดเริ่มต้นของ 99 และสิ้นสุด 115 นี่เป็นเพียงเพื่อความสอดคล้องกับไฟล์ Bigwig พื้นฐานและอาจเปลี่ยนแปลงในอนาคต
Pybigwig ยังมีอยู่ในแพ็คเกจใน Galaxy คุณสามารถค้นหาได้ในเครื่องมือและ IUC กำลังโฮสต์นิยาม XML ของสิ่งนี้ใน GitHub
