השתמש ב- picard FilterSamReads http://broadinstitute.github.io/picard/command-line-overview.html#FilterSamReads

READ_LIST_FILE (קובץ) רשימת קריאה קובץ המכיל קריאות שייכלל או לא ייכלל בקובץ OUTPUT SAM או BAM. ערך ברירת מחדל: null.

(1) זה לא יהיה מהיר במיוחד, אבל אתה יכול לספק ל- grep קובץ QNAMES.

  samtools view file.bam | grep -f 'qnames.txt > subset.sam  

where qnames.txt has

  דוגמא: QNAME1EXAMPLE: QNAME2EXAMPLE: QNAME3EXAMPLE: QNAME4EXAMPLE: QNAME5  קוד> 

(2) זה יהיה קצת יותר מסובך, אך האם אתה יכול לתת דוגמה שבה ל- grep יכול להיות QNAME הנכון?

(3) כדי לשמור על כותרת BAM הייתי משתמש ב samtools -H file.bam > header.txt כדי להשיג את הכותרת ואז חתול את קובץ הכותרת עם הקובץ grep'ed

להלן קטע קטן של קוד פיתון שמראה דרך אחת נאיבית אך ניתנת לניהול (הערה, אני מצפה ש- grep יהיה מהיר יותר, אם כי לגרום לו להוציא את הכותרת יהיה מעצבן):

  ייבוא ​​pysamqnames = set (...) # שמות קריאה עבור לכאןbam = pysam.AlignmentFile ("איזה קובץ קלט.באם") obam = pysam. AlignmentFile ("קובץ פלט כלשהו. Bam", "w", תבנית = bam) עבור b ב bam.fetch (עד_ eof = נכון): אם b.query_name בשמות שמות: obam.write (b) bam.close () obam. סגור ()  

עכשיו זה יעבוד, אבל אם יש לך הרבה שמות לעיין אז זה יהיה די איטי (כמו ש- grep -f. .. ). אם עליך לבחור רשימה גדולה של קריאות לפי שם, אז אסטרטגיה יעילה יותר היא:

1. שם מיון קובץ ה- BAM ( samtools sort -n או picard) , שניתן לעשות עם מספר שרשורים.
2. מיין את רשימת שמות הקריאה שיתאימו (זה יהיה מסובך יותר ממה שאפשר היה לחשוב בתמימות, מכיוון ש- samtools ו- picard יתנו שם אחרת, אז הקפד לשים חשבו על זה).
3. בצע את אותו סוג של איטרציה כמו לעיל, אך השווה רק לשם הקריאה העליון ברשימה שלך משלב 2 (הסרת אותו אלמנט מהערימה לאחר התאמה).

ביצוע grep ב $ n $ יישור ל $ m $ שמות יתן לך פעולות $ O (mn) $ . אם תחליט למיין תחילה את שמות ה- BAM והשמות שלך, תוכל לצמצם אותם ל $ O (\ min \ {m, n \}) $ על ידי הקפצת קריאות ושמות qn. ערימותיהם. אתה גם לא קורא את כל ה- BAM בבת אחת ככה, מה שמגביל את צריכת הזיכרון. אני די בטוח שגם פיקארד עושה משהו כזה.

הנה דוגמה ב- Python שכתבתי בדיוק למטרה זו:

ניתן להשתמש ב SAMsift:

  samsift -i file.bam -0 'q = open ("qnames .txt "). קרא (). splitlines () '-f' QNAME ב- q ' 

-i file.bam מציין את קובץ ה- BAM הקלט. -0 'q = open ("qnames.txt"). קרא (). splitlines ()' טוען את רשימת ה- QNAME שלך. -f 'QNAME in q' מציין את המסנן ליישור.

עם SAMsift אפשר להשתמש בכל ביטוי Python לסינון. החיסרון הוא שהוא איטי יותר מהכלים האחרים.

מרבית הפתרונות שפורסמו עד כה הם איטיים או יכולים להניב תוצאות מזויפות במקרים מסוימים. (לא בדקתי את פיקארד, מה שאני מניח שהוא יעבוד כמתוכנן. אבל באופן אישי אני נוטה לרתוע כשאצטרך להקים מכשיר JVM!)

בוא ניגש לבעיה מזווית אחרת:

לפונקציות hash טובות יש מורכבות ~ O (1) לבדיקה. אם אנו יכולים ליצור טבלת hash של ה- QNAME שלנו, אז חיפוש קובץ BAM צריך להיות O (n) במספר הרשומות. מפתחי ה- Python גאים מאוד בפונקציית ה- hash ב- cpython, אז בואו נשתמש בזה:

  #! / Usr / bin / env python3import syswith open (sys.argv [1], 'r' ) כקובץ אינדקס: ids = set (l.rstrip ('\ r \ n') עבור l בקובץ index) עבור שורה ב- sys.stdin: qname, _ = line.split ('\ t', 1) אם qname בתעודות : sys.stdout.write (line)  

עם ערכת בדיקה המורכבת מקובץ BAM ~ 3GB ו qnames.txt עם ~ 65,000 ערכים, פועל samtools להציג input.bam | ./idfilter.py qnames.txt לוקח בערך 6.5 שניות בשרת הבדיקה שלי. לשם השוואה, צנרת ה- SAM ל fgrep -f qnames.txt (באמצעות GNU grep) אורכת כ- 8.5 שניות (אך עשויה להניב תוצאות מזויפות).

(מדוע אולי grep מניב תוצאות מזויפות? אם יש לך foo QNAME ב- qnames.txt , אך קובץ ה- BAM שלך מכיל גם fooX ו- fooY , כולם יותאמו על ידי fgrep . הפיתרון הוא לעגן כל אחד מדפוסי החיפוש שלך בין תחילת השורה לכרטיסייה, למשל ^ foo \ t , אך אז עליך להשתמש ב grep סטנדרטי אשר יצטרך לבנות NFA עבור כל תבנית, ולא ב fgrep אשר מיישם את Aho-Corasick, ו- תהיה סדרי גודל איטיים יותר.)

כתבתי מחדש את סקריפט הפיתון הקטן שלי ב- Rust באמצעות std :: אוספים :: HashMap לפונקציית ה- hash, ובאמצעות הארגז rust_htslib לקרוא ולכתוב היישר מ- BAM, וזה עוד מהר פי שניים מהפייתון, גם כשקוראים וכותבים BAM. עם זאת, החלודה שלי לא ממש מתאימה לצריכה ציבורית ...