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

"רק תראה לי את הפתרון!" אם אתה רק רוצה לראות מאיפה באתי, אז דלג לפרק "להיות יותר יצירתי", אבל אני חושב שזה מעניין ושימושי יותר לקרוא על כישלון.

לאחרונה הוטל עלי להקים תהליך לעיבוד נפח גדול של רצפי DNA גולמיים (טכנית שבב SNP). הצורך היה להשיג במהירות נתונים על מיקום גנטי נתון (הנקרא SNP) עבור מודלים שלאחר מכן ומשימות אחרות. באמצעות R ו-AWK, הצלחתי לנקות ולארגן נתונים בצורה טבעית, מה שהאיץ מאוד את עיבוד השאילתות. זה לא היה קל עבורי ודרש איטרציות רבות. מאמר זה יעזור לך להימנע מכמה מהטעויות שלי ויראה לך מה הגעתי אליו.

ראשית, כמה הסברים מבוא.

נתונים

המרכז לעיבוד מידע גנטי באוניברסיטה סיפק לנו נתונים בצורה של TSV של 25 TB. קיבלתי אותם מחולקים ל-5 חבילות, דחוסות על ידי Gzip, שכל אחת מהן הכילה כ-240 קבצים של ארבעה גיגה. כל שורה הכילה נתונים עבור SNP אחד מאדם אחד. בסך הכל, הועברו נתונים על ~2,5 מיליון SNP ו~60 אלף אנשים. בנוסף למידע SNP, הקבצים הכילו עמודות רבות עם מספרים המשקפים מאפיינים שונים, כגון עוצמת קריאה, תדירות של אללים שונים וכו'. בסך הכל היו כ-30 עמודות עם ערכים ייחודיים.

יעד

כמו בכל פרויקט ניהול נתונים, הדבר החשוב ביותר היה לקבוע כיצד ייעשה שימוש בנתונים. במקרה הזה אנו נבחר בעיקר מודלים וזרימות עבודה עבור SNP על בסיס SNP. כלומר, נצטרך רק נתונים על SNP אחד בכל פעם. הייתי צריך ללמוד איך לאחזר את כל הרשומות הקשורות לאחד מ-2,5 מיליון SNPs בקלות, במהירות ובזול ככל האפשר.

איך לא לעשות את זה

אם לצטט קלישאה מתאימה:

לא נכשלתי אלף פעמים, פשוט גיליתי אלף דרכים להימנע מניתוח של חבורה של נתונים בפורמט ידידותי לשאילתות.

ניסיון ראשון

מה למדתי: אין דרך זולה לנתח 25 TB בכל פעם.

לאחר שלמדתי את הקורס "שיטות מתקדמות לעיבוד ביג דאטה" באוניברסיטת ונדרבילט, הייתי בטוח שהטריק נמצא בקופה. זה ייקח כנראה שעה או שעתיים להגדיר את שרת Hive לרוץ על כל הנתונים ולדווח על התוצאה. מכיוון שהנתונים שלנו מאוחסנים ב-AWS S3, השתמשתי בשירות אתנה, המאפשר לך להחיל שאילתות Hive SQL על נתוני S3. אינך צריך להקים/להעלות אשכול Hive, ואתה גם משלם רק עבור הנתונים שאתה מחפש.

לאחר שהראיתי לאתנה את הנתונים שלי ואת הפורמט שלהם, הרצתי כמה בדיקות עם שאילתות כמו זה:

select * from intensityData limit 10;

וקיבל במהירות תוצאות מובנות היטב. מוּכָן.

עד שניסינו להשתמש בנתונים בעבודתנו...

התבקשתי לשלוף את כל מידע SNP כדי לבדוק את הדגם. הרצתי את השאילתה:

select * from intensityData 
where snp = 'rs123456';

...והחל לחכות. לאחר שמונה דקות ויותר מ-4 TB של נתונים מבוקשים, קיבלתי את התוצאה. אתנה גובה לפי נפח הנתונים שנמצאו, $5 לטרה-בייט. אז הבקשה הבודדת הזו עלתה 20 דולר ושמונה דקות של המתנה. כדי להפעיל את המודל על כל הנתונים, היינו צריכים לחכות 38 שנים ולשלם 50 מיליון דולר, ברור שזה לא התאים לנו.

היה צורך להשתמש בפרקט...

מה למדתי: היזהר עם גודל קבצי הפרקט והארגון שלהם.

תחילה ניסיתי לתקן את המצב על ידי המרת כל TSVs ל קבצי פרקט. הם נוחים לעבודה עם מערכי נתונים גדולים מכיוון שהמידע בהם מאוחסן בצורה עמודת: כל עמודה נמצאת בקטע זיכרון/דיסק משלה, בניגוד לקבצי טקסט, שבהם שורות מכילות אלמנטים של כל עמודה. ואם אתה צריך למצוא משהו, אז פשוט קרא את הטור הנדרש. בנוסף, כל קובץ מאחסן טווח של ערכים בעמודה, כך שאם הערך שאתה מחפש אינו בטווח העמודה, Spark לא יבזבז זמן בסריקת הקובץ כולו.

רצתי משימה פשוטה דבק AWS להמיר את ה-TSVs שלנו ל-Parquet והורדת הקבצים החדשים לתוך Athena. זה לקח בערך 5 שעות. אבל כשהרצתי את הבקשה, זה לקח בערך אותו פרק זמן וקצת פחות כסף להשלים. העובדה היא שספארק, בניסיון לייעל את המשימה, פשוט פירק נתח TSV אחד ושם אותו בנתח פרקט משלו. ומכיוון שכל נתח היה גדול מספיק כדי להכיל את כל הרשומות של אנשים רבים, כל קובץ הכיל את כל ה-SNPs, אז Spark היה צריך לפתוח את כל הקבצים כדי לחלץ את המידע שהוא צריך.

מעניין שסוג הדחיסה המומלץ (והמומלץ) של Parquet, Snappy, אינו ניתן לפיצול. לכן, כל מבצע היה תקוע במשימה של פירוק והורדת מערך הנתונים המלא של 3,5 GB.

בואו נבין את הבעיה

מה למדתי: המיון קשה, במיוחד אם הנתונים מופצים.

נראה לי שעכשיו הבנתי את מהות הבעיה. הייתי צריך רק למיין את הנתונים לפי עמודת SNP, לא לפי אנשים. לאחר מכן, מספר SNPs יאוחסנו בנתח נתונים נפרד, ואז הפונקציה "חכמה" של Parquet "פתוח רק אם הערך נמצא בטווח" תראה את עצמה במלוא תפארתה. לרוע המזל, מיון של מיליארדי שורות הפזורות על פני אשכול התברר כמשימה קשה.

אני לוקח שיעור אלגוריתמים בקולג': "אוף, לאף אחד לא אכפת מהמורכבות החישובית של כל אלגוריתמי המיון האלה"

אני מנסה למיין על עמודה ב-20TB #לְעוֹרֵר טבלה: "למה זה לוקח כל כך הרבה זמן?" # DataScience מאבקים.

— ניק סטרייר (@NicholasStrayer) במרץ 11, 2019

AWS בהחלט לא רוצה להנפיק החזר בגלל הסיבה "אני תלמיד מוסח". אחרי שהרצתי מיון ב- Amazon Glue, הוא רץ במשך יומיים וקרס.

מה לגבי חלוקה למחיצות?

מה למדתי: מחיצות ב-Spark חייבות להיות מאוזנות.

ואז הגעתי לרעיון של חלוקת נתונים בכרומוזומים. יש 23 מהם (ועוד כמה אם לוקחים בחשבון DNA מיטוכונדריאלי ואזורים לא ממופים).
זה יאפשר לך לפצל את הנתונים לנתחים קטנים יותר. אם תוסיף רק שורה אחת לפונקציית הייצוא Spark בסקריפט Glue partition_by = "chr", אז יש לחלק את הנתונים לדליים.

הגנום מורכב משברים רבים הנקראים כרומוזומים.

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

מה אם נחלק אותו לחלקים קטנים עוד יותר?

מה למדתי: לעולם אל תנסה לעשות 2,5 מיליון מחיצות בכלל.

החלטתי ללכת עד הסוף ולחלק כל SNP. זה הבטיח שהמחיצות יהיו בגודל שווה. זה היה רעיון גרוע. השתמשתי בדבק והוספתי קו תמים partition_by = 'snp'. המשימה התחילה והחלה לצאת לפועל. יום לאחר מכן בדקתי וראיתי שעדיין לא כתוב כלום ל-S3, אז הרגתי את המשימה. זה נראה כאילו Glue כתב קבצי ביניים למיקום נסתר ב-S3, הרבה קבצים, אולי כמה מיליון. כתוצאה מכך, הטעות שלי עלתה יותר מאלף דולר ולא מצאה חן בעיני המנטור שלי.

חלוקה + מיון

מה למדתי: המיון עדיין קשה, וכך גם כוונון Spark.

הניסיון האחרון שלי לחלוקה כלל חלוקה של הכרומוזומים ולאחר מכן מיון כל מחיצה. בתיאוריה, זה יזרז כל שאילתה מכיוון שנתוני ה-SNP הרצויים היו צריכים להיות בתוך כמה נתחי פרקט בטווח נתון. לרוע המזל, מיון אפילו נתונים מחולקים התברר כמשימה קשה. כתוצאה מכך, עברתי ל-EMR עבור אשכול מותאם אישית והשתמשתי בשמונה מופעים רבי עוצמה (C5.4xl) וב-Sparklyr כדי ליצור זרימת עבודה גמישה יותר...

# Sparklyr snippet to partition by chr and sort w/in partition
# Join the raw data with the snp bins
raw_data
  group_by(chr) %>%
  arrange(Position) %>% 
  Spark_write_Parquet(
    path = DUMP_LOC,
    mode = 'overwrite',
    partition_by = c('chr')
  )

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

עדכון: אז זה מתחיל. pic.twitter.com/agY4GU2ru5

— ניק סטרייר (@NicholasStrayer) מאי 15, 2019

אני נהיה יותר יצירתי

מה למדתי: לפעמים נתונים מיוחדים דורשים פתרונות מיוחדים.

לכל SNP יש ערך מיקום. זהו מספר המתאים למספר הבסיסים לאורך הכרומוזום שלו. זוהי דרך נחמדה וטבעית לארגן את הנתונים שלנו. בהתחלה רציתי לחלק לפי אזורים של כל כרומוזום. לדוגמה, עמדות 1 - 2000, 2001 - 4000 וכו'. אבל הבעיה היא ש-SNPs אינם מפוזרים באופן שווה על פני הכרומוזומים, ולכן גדלי הקבוצה ישתנו מאוד.

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

בסופו של דבר, הכנתי קבוצות (פח) של 75 SNP, הסיבה תוסבר בהמשך.

snp_to_bin <- unique_snps %>% 
  group_by(chr) %>% 
  arrange(position) %>% 
  mutate(
    rank = 1:n()
    bin = floor(rank/snps_per_bin)
  ) %>% 
  ungroup()

נסה קודם עם ספארק

מה למדתי: צבירת הניצוץ היא מהירה, אבל החלוקה עדיין יקרה.

רציתי לקרוא את מסגרת הנתונים הקטנה הזו (2,5 מיליון שורות) לתוך Spark, לשלב אותה עם הנתונים הגולמיים, ואז לחלק אותה לפי העמודה החדשה שנוספה bin.

# Join the raw data with the snp bins
data_w_bin <- raw_data %>%
  left_join(sdf_broadcast(snp_to_bin), by ='snp_name') %>%
  group_by(chr_bin) %>%
  arrange(Position) %>% 
  Spark_write_Parquet(
    path = DUMP_LOC,
    mode = 'overwrite',
    partition_by = c('chr_bin')
  )

השתמשתי sdf_broadcast(), אז Spark יודע שהוא צריך לשלוח את מסגרת הנתונים לכל הצמתים. זה שימושי אם הנתונים קטנים בגודלם ונדרשים עבור כל המשימות. אחרת, Spark מנסה להיות חכם ומפיץ נתונים לפי הצורך, מה שעלול לגרום להאטות.

ושוב, הרעיון שלי לא עבד: המשימות עבדו במשך זמן מה, השלימו את האיחוד, ואז, כמו המבצעים שהושקו על ידי חלוקה, הם החלו להיכשל.

הוספת AWK

מה למדתי: אל תישן כאשר מלמדים אותך את היסודות. בטח מישהו כבר פתר את הבעיה שלך בשנות ה-1980.

עד לנקודה זו, הסיבה לכל הכישלונות שלי עם Spark הייתה ערבוביית הנתונים באשכול. אולי ניתן לשפר את המצב עם טיפול מקדים. החלטתי לנסות לפצל את נתוני הטקסט הגולמיים לעמודות של כרומוזומים, אז קיוויתי לספק ל-Spark נתונים "מחולקים מראש".

חיפשתי ב-StackOverflow כיצד לפצל לפי ערכי עמודות ומצאתי תשובה כל כך נהדרת. עם AWK אתה יכול לפצל קובץ טקסט לפי ערכי עמודות על ידי כתיבתו בסקריפט במקום שליחת התוצאות אל stdout.

כתבתי תסריט של Bash כדי לנסות אותו. הורד את אחד מה-TSVs הארוזים, ואז פרק אותו באמצעות gzip ונשלח ל awk.

gzip -dc path/to/chunk/file.gz |
awk -F 't' 
'{print $1",..."$30">"chunked/"$chr"_chr"$15".csv"}'

זה עבד!

מילוי הליבות

מה למדתי: gnu parallel - זה דבר קסום, כולם צריכים להשתמש בו.

ההפרדה הייתה די איטית וכשהתחלתי htopכדי לבדוק את השימוש במופע חזק (ויקר) של EC2, התברר שאני משתמש רק בליבה אחת ובערך 200 מגה-בייט של זיכרון. כדי לפתור את הבעיה ולא להפסיד הרבה כסף, היינו צריכים להבין איך מקבילים את העבודה. למרבה המזל, בספר מדהים לחלוטין מדע נתונים בשורת הפקודה מצאתי פרק מאת ג'רון יאנסנס על הקבלה. ממנו למדתי על gnu parallel, שיטה גמישה מאוד להטמעת ריבוי פתילים ביוניקס.

כשהתחלתי את החלוקה באמצעות התהליך החדש, הכל היה בסדר, אבל עדיין היה צוואר בקבוק - הורדת אובייקטי S3 לדיסק לא הייתה מהירה במיוחד ולא הייתה מקבילה לחלוטין. כדי לתקן את זה, עשיתי את זה:

1. גיליתי שאפשר ליישם את שלב ההורדה של S3 ישירות בצנרת, ולבטל לחלוטין את אחסון הביניים בדיסק. זה אומר שאני יכול להימנע מכתיבת נתונים גולמיים לדיסק ולהשתמש אפילו באחסון קטן יותר, ולכן זול יותר, ב-AWS.
2. קְבוּצָה aws configure set default.s3.max_concurrent_requests 50 הגדיל מאוד את מספר השרשורים שבהם משתמש AWS CLI (כברירת מחדל יש 10).
3. עברתי למופע EC2 מותאם למהירות רשת, עם האות n בשם. גיליתי שאובדן כוח העיבוד בעת שימוש ב-n-מופעים יותר מפיצוי על ידי העלייה במהירות הטעינה. לרוב המשימות השתמשתי ב-c5n.4xl.
4. השתנה gzip על pigz, זהו כלי gzip שיכול לעשות דברים מגניבים כדי להקביל את המשימה הלא-מקבילה בתחילה של פירוק קבצים (זה הכי פחות עזר).

# Let S3 use as many threads as it wants
aws configure set default.s3.max_concurrent_requests 50

for chunk_file in $(aws s3 ls $DATA_LOC | awk '{print $4}' | grep 'chr'$DESIRED_CHR'.csv') ; do

        aws s3 cp s3://$batch_loc$chunk_file - |
        pigz -dc |
        parallel --block 100M --pipe  
        "awk -F 't' '{print $1",..."$30">"chunked/{#}_chr"$15".csv"}'"

       # Combine all the parallel process chunks to single files
        ls chunked/ |
        cut -d '_' -f 2 |
        sort -u |
        parallel 'cat chunked/*_{} | sort -k5 -n -S 80% -t, | aws s3 cp - '$s3_dest'/batch_'$batch_num'_{}'
        
         # Clean up intermediate data
       rm chunked/*
done

השלבים הללו משולבים זה בזה כדי שהכל יעבוד מהר מאוד. על ידי הגדלת מהירויות ההורדה וביטול כתיבת דיסקים, יכולתי כעת לעבד חבילה של 5 טרה-בייט תוך מספר שעות בלבד.

אין דבר מתוק יותר מלראות את כל הליבות שאתה משלם עליהן ב-AWS בשימוש. הודות ל-gnu-parallel אני יכול לפתוח ולפצל CSV של 19 גיגה בדיוק באותה מהירות שאני יכול להוריד אותו. אפילו לא הצלחתי להשיג ניצוץ כדי להפעיל את זה. # DataScience # לינוקס pic.twitter.com/Nqyba2zqEk

— ניק סטרייר (@NicholasStrayer) מאי 17, 2019

הציוץ הזה היה צריך להזכיר את 'TSV'. אבוי.

שימוש בנתונים שמנותחו לאחרונה

מה למדתי: Spark אוהב נתונים לא דחוסים ולא אוהב שילוב מחיצות.

כעת הנתונים היו ב-S3 בפורמט לא ארוז (לקרא: משותף) ומסודר למחצה, ויכולתי לחזור שוב ל-Spark. חיכתה לי הפתעה: שוב לא הצלחתי להשיג את מה שרציתי! היה קשה מאוד לומר ל-Spark בדיוק איך הנתונים מחולקים. וגם כשעשיתי את זה, התברר שיש יותר מדי מחיצות (95 אלף), וכשהשתמשתי coalesce צמצם את מספרם לגבולות סבירים, זה הרס את החלוקה שלי. אני בטוח שאפשר לתקן את זה, אבל אחרי כמה ימים של חיפושים לא הצלחתי למצוא פתרון. בסופו של דבר סיימתי את כל המשימות ב-Spark, למרות שזה לקח זמן וקבצי ה-Parquet המפוצלים שלי לא היו קטנים במיוחד (~200 KB). עם זאת, הנתונים היו היכן שהיה צורך.

קטן מדי ולא אחיד, נפלא!

בדיקת שאילתות Spark מקומיות

מה למדתי: ל-Spark יש יותר מדי תקורה בעת פתרון בעיות פשוטות.

על ידי הורדת הנתונים בפורמט חכם, הצלחתי לבדוק את המהירות. הגדר סקריפט R להפעלת שרת Spark מקומי, ולאחר מכן טען מסגרת נתונים של Spark מהאחסון הקבוצתי של Parquet (סל) שצוין. ניסיתי לטעון את כל הנתונים אך לא הצלחתי לגרום ל-Sparklyr לזהות את החלוקה.

sc <- Spark_connect(master = "local")

desired_snp <- 'rs34771739'

# Start a timer
start_time <- Sys.time()

# Load the desired bin into Spark
intensity_data <- sc %>% 
  Spark_read_Parquet(
    name = 'intensity_data', 
    path = get_snp_location(desired_snp),
    memory = FALSE )

# Subset bin to snp and then collect to local
test_subset <- intensity_data %>% 
  filter(SNP_Name == desired_snp) %>% 
  collect()

print(Sys.time() - start_time)

ההוצאה להורג ארכה 29,415 שניות. הרבה יותר טוב, אבל לא טוב מדי לבדיקה המונית של כל דבר. בנוסף, לא יכולתי להאיץ את העניינים עם שמירה במטמון, כי כשניסיתי לאחסן מסגרת נתונים בזיכרון, Spark תמיד קרס, גם כשהקצאתי יותר מ-50 GB של זיכרון למערך נתונים ששקל פחות מ-15.

חזור ל-AWK

מה למדתי: מערכים אסוציאטיביים ב-AWK יעילים מאוד.

הבנתי שאני יכול להשיג מהירויות גבוהות יותר. נזכרתי בזה בצורה נפלאה מדריך AWK מאת ברוס בארנט קראתי על תכונה מגניבה בשם "מערכים אסוציאטיביים" בעיקרו של דבר, אלו הם צמדי מפתח-ערך, שמשום מה נקראו אחרת ב-AWK, ולכן איכשהו לא חשבתי עליהם הרבה. רומן צ'פליאקה נזכר שהמונח "מערכים אסוציאטיביים" עתיק בהרבה מהמונח "זוג מפתח-ערך". גם אם אתה חפש את ערך המפתח ב-Google Ngram, לא תראה שם את המונח הזה, אבל תמצא מערכים אסוציאטיביים! בנוסף, "צמד ערך-מפתח" משויך לרוב לבסיסי נתונים, ולכן הרבה יותר הגיוני להשוות אותו עם hashmap. הבנתי שאני יכול להשתמש במערכים האסוציאטיביים האלה כדי לשייך את ה-SNPs שלי לטבלת bin ונתונים גולמיים מבלי להשתמש ב-Spark.

לשם כך, בסקריפט AWK השתמשתי בבלוק BEGIN. זהו קטע קוד שמבוצע לפני ששורת הנתונים הראשונה מועברת לגוף הראשי של הסקריפט.

join_data.awk
BEGIN {
  FS=",";
  batch_num=substr(chunk,7,1);
  chunk_id=substr(chunk,15,2);
  while(getline < "snp_to_bin.csv") {bin[$1] = $2}
}
{
  print $0 > "chunked/chr_"chr"_bin_"bin[$1]"_"batch_num"_"chunk_id".csv"
}

קבוצה while(getline...) טען את כל השורות מקבוצת ה-CSV (bin), הגדר את העמודה הראשונה (שם SNP) כמפתח עבור המערך האסוציאטיבי bin והערך השני (קבוצה) כערך. ואז בבלוק { }, שמתבצע בכל שורות הקובץ הראשי, כל שורה נשלחת לקובץ הפלט, שמקבל שם ייחודי בהתאם לקבוצה שלו (bin): ..._bin_"bin[$1]"_....

משתנים batch_num и chunk_id תאמו את הנתונים שסיפקו הצינור, תוך הימנעות ממצב מירוץ, וכל חוט ביצוע פועל parallel, כתב לקובץ ייחודי משלו.

מכיוון שפיזרתי את כל הנתונים הגולמיים לתיקיות על כרומוזומים שנשארו מהניסוי הקודם שלי עם AWK, עכשיו יכולתי לכתוב סקריפט Bash אחר כדי לעבד כרומוזום אחד בכל פעם ולשלוח נתונים מחולקים עמוקים יותר ל-S3.

DESIRED_CHR='13'

# Download chromosome data from s3 and split into bins
aws s3 ls $DATA_LOC |
awk '{print $4}' |
grep 'chr'$DESIRED_CHR'.csv' |
parallel "echo 'reading {}'; aws s3 cp "$DATA_LOC"{} - | awk -v chr=""$DESIRED_CHR"" -v chunk="{}" -f split_on_chr_bin.awk"

# Combine all the parallel process chunks to single files and upload to rds using R
ls chunked/ |
cut -d '_' -f 4 |
sort -u |
parallel "echo 'zipping bin {}'; cat chunked/*_bin_{}_*.csv | ./upload_as_rds.R '$S3_DEST'/chr_'$DESIRED_CHR'_bin_{}.rds"
rm chunked/*

לתסריט יש שני חלקים parallel.

בסעיף הראשון קוראים נתונים מכל הקבצים המכילים מידע על הכרומוזום הרצוי, לאחר מכן נתונים אלו מופצים על פני שרשורים, אשר מפיצים את הקבצים לקבוצות המתאימות (bin). כדי להימנע מתנאי מרוץ כאשר שרשורים מרובים כותבים לאותו קובץ, AWK מעביר את שמות הקבצים כדי לכתוב נתונים למקומות שונים, למשל. chr_10_bin_52_batch_2_aa.csv. כתוצאה מכך נוצרים קבצים קטנים רבים בדיסק (בשביל זה השתמשתי בנפחי EBS של טרה-בייט).

מסוע מהחלק השני parallel עובר על הקבוצות (bin) ומשלב את הקבצים האישיים שלהם ל-CSV משותף c catולאחר מכן שולח אותם לייצוא.

שידור ב-R?

מה למדתי: אתה יכול ליצור קשר stdin и stdout מסקריפט R, ולכן משתמשים בו בצנרת.

אולי שמת לב לשורה הזו בתסריט הבאש שלך: ...cat chunked/*_bin_{}_*.csv | ./upload_as_rds.R.... זה מתרגם את כל קבצי הקבוצה המשורשרים (bin) לסקריפט R למטה. {} היא טכניקה מיוחדת parallel, אשר מכניס כל נתונים שהוא שולח לזרם שצוין ישירות לתוך הפקודה עצמה. אוֹפְּצִיָה {#} מספק מזהה שרשור ייחודי, ו {%} מייצג את מספר משבצת העבודה (חוזר, אך לעולם לא בו-זמנית). רשימה של כל האפשרויות ניתן למצוא ב תיעוד.

#!/usr/bin/env Rscript
library(readr)
library(aws.s3)

# Read first command line argument
data_destination <- commandArgs(trailingOnly = TRUE)[1]

data_cols <- list(SNP_Name = 'c', ...)

s3saveRDS(
  read_csv(
        file("stdin"), 
        col_names = names(data_cols),
        col_types = data_cols 
    ),
  object = data_destination
)

כאשר משתנה file("stdin") מועבר ל readr::read_csv, הנתונים המתורגמים לסקריפט R נטען לתוך מסגרת, אשר לאחר מכן בטופס .rds-קובץ באמצעות aws.s3 נכתב ישירות ל-S3.

RDS הוא משהו כמו גרסה זוטרת של פרקט, ללא התלבטויות של אחסון רמקולים.

לאחר סיום התסריט של Bash קיבלתי חבילה .rds-קבצים הממוקמים ב-S3, מה שאפשר לי להשתמש בדחיסה יעילה ובסוגים מובנים.

למרות השימוש בבלם R, הכל עבד מהר מאוד. באופן לא מפתיע, החלקים של R שקוראים וכותבים נתונים מותאמים מאוד. לאחר בדיקה על כרומוזום אחד בגודל בינוני, העבודה הושלמה על מופע C5n.4xl תוך כשעתיים.

מגבלות S3

מה למדתי: הודות ליישום נתיב חכם, S3 יכול להתמודד עם קבצים רבים.

דאגתי אם S3 יוכל לטפל בקבצים הרבים שהועברו אליו. אני יכול לגרום לשמות הקבצים להיות הגיוניים, אבל איך S3 יחפש אותם?

תיקיות ב-S3 הן רק לראווה, למעשה המערכת לא מעוניינת בסמל /. מדף השאלות הנפוצות של S3.

נראה ש-S3 מייצג את הנתיב לקובץ מסוים כמפתח פשוט במעין טבלת hash או מסד נתונים מבוסס מסמכים. ניתן להתייחס לדלי כעל טבלה, וקבצים יכולים להיחשב כרשומות בטבלה זו.

מכיוון שמהירות ויעילות חשובות לרווח באמזון, אין זה מפתיע שמערכת המפתח כקובץ נתיב זו אופטימלית בצורה מטורפת. ניסיתי למצוא איזון: כדי שלא אצטרך לבצע הרבה בקשות get, אלא שהבקשות בוצעו במהירות. התברר שעדיף להכין כ-20 אלף קבצי bin. אני חושב שאם נמשיך לבצע אופטימיזציה, נוכל להשיג עלייה במהירות (לדוגמה, יצירת דלי מיוחד רק עבור נתונים, ובכך להקטין את גודל טבלת החיפוש). אבל לא היה זמן או כסף לניסויים נוספים.

מה לגבי תאימות צולבת?

מה למדתי: הסיבה מספר אחת לבזבוז זמן היא אופטימיזציה של שיטת האחסון שלך בטרם עת.

בשלב זה, חשוב מאוד לשאול את עצמך: "מדוע להשתמש בפורמט קובץ קנייני?" הסיבה נעוצה במהירות הטעינה (לטעינת קבצי CSV ב-gzip לקח פי 7 יותר זמן) ובתאימות עם זרימות העבודה שלנו. אני עשוי לשקול מחדש אם R יכול לטעון בקלות קבצי פרקט (או חץ) ללא עומס Spark. כולם במעבדה שלנו משתמשים ב-R, ואם אני צריך להמיר את הנתונים לפורמט אחר, עדיין יש לי את נתוני הטקסט המקוריים, אז אני יכול פשוט להפעיל את הצינור שוב.

חלוקת עבודה

מה למדתי: אל תנסה לייעל עבודות באופן ידני, תן ​​למחשב לעשות זאת.

ניפוי באגים בזרימת העבודה בכרומוזום אחד, עכשיו אני צריך לעבד את כל הנתונים האחרים.
רציתי להעלות כמה מופעים של EC2 להמרה, אבל במקביל פחדתי לקבל עומס מאוד לא מאוזן על פני עבודות עיבוד שונות (בדיוק כפי ש-Spark סבל ממחיצות לא מאוזנות). בנוסף, לא הייתי מעוניין להעלות מופע אחד לכל כרומוזום, כי לחשבונות AWS יש מגבלת ברירת מחדל של 10 מופעים.

ואז החלטתי לכתוב סקריפט ב-R כדי לייעל עבודות עיבוד.

ראשית, ביקשתי מ-S3 לחשב כמה שטח אחסון תופס כל כרומוזום.

library(aws.s3)
library(tidyverse)

chr_sizes <- get_bucket_df(
  bucket = '...', prefix = '...', max = Inf
) %>% 
  mutate(Size = as.numeric(Size)) %>% 
  filter(Size != 0) %>% 
  mutate(
    # Extract chromosome from the file name 
    chr = str_extract(Key, 'chr.{1,4}.csv') %>%
             str_remove_all('chr|.csv')
  ) %>% 
  group_by(chr) %>% 
  summarise(total_size = sum(Size)/1e+9) # Divide to get value in GB



# A tibble: 27 x 2
   chr   total_size
   <chr>      <dbl>
 1 0           163.
 2 1           967.
 3 10          541.
 4 11          611.
 5 12          542.
 6 13          364.
 7 14          375.
 8 15          372.
 9 16          434.
10 17          443.
# … with 17 more rows

אחר כך כתבתי פונקציה שלוקחת את הגודל הכולל, מערבבת את סדר הכרומוזומים, מחלקת אותם לקבוצות num_jobs ומספר לך עד כמה הגדלים שונים של כל עבודות העיבוד.

num_jobs <- 7
# How big would each job be if perfectly split?
job_size <- sum(chr_sizes$total_size)/7

shuffle_job <- function(i){
  chr_sizes %>%
    sample_frac() %>% 
    mutate(
      cum_size = cumsum(total_size),
      job_num = ceiling(cum_size/job_size)
    ) %>% 
    group_by(job_num) %>% 
    summarise(
      job_chrs = paste(chr, collapse = ','),
      total_job_size = sum(total_size)
    ) %>% 
    mutate(sd = sd(total_job_size)) %>% 
    nest(-sd)
}

shuffle_job(1)



# A tibble: 1 x 2
     sd data            
  <dbl> <list>          
1  153. <tibble [7 × 3]>

אחר כך עברתי אלף דשדושים באמצעות purrr ובחרתי את הטוב ביותר.

1:1000 %>% 
  map_df(shuffle_job) %>% 
  filter(sd == min(sd)) %>% 
  pull(data) %>% 
  pluck(1)

אז בסופו של דבר קיבלתי קבוצה של משימות שהיו מאוד דומות בגודלן. ואז כל מה שנותר היה לעטוף את התסריט הקודם שלי בשש בלופ גדול for. כתיבת האופטימיזציה הזו ארכה כ-10 דקות. וזה הרבה פחות ממה שהייתי מוציא על יצירה ידנית של משימות אם הן לא היו מאוזנות. לכן, אני חושב שצדקתי עם האופטימיזציה הראשונית הזו.

for DESIRED_CHR in "16" "9" "7" "21" "MT"
do
# Code for processing a single chromosome
fi

בסוף אני מוסיף את פקודת הכיבוי:

sudo shutdown -h now

... והכל הסתדר! באמצעות ה-AWS CLI, העליתי מופעים באמצעות האפשרות user_data נתן להם תסריטים של Bash של המשימות שלהם לעיבוד. הם רצו ונכבו אוטומטית, אז לא שילמתי עבור כוח עיבוד נוסף.

aws ec2 run-instances ...
--tag-specifications "ResourceType=instance,Tags=[{Key=Name,Value=<<job_name>>}]" 
--user-data file://<<job_script_loc>>

בואו לארוז!

מה למדתי: ה-API צריך להיות פשוט למען הקלות והגמישות בשימוש.

לבסוף קיבלתי את הנתונים במקום ובצורה הנכונים. כל שנותר היה לפשט את תהליך השימוש בנתונים ככל האפשר כדי להקל על הקולגות שלי. רציתי ליצור API פשוט ליצירת בקשות. אם בעתיד אחליט לעבור מ .rds לקבצי פרקט, אז זו צריכה להיות בעיה עבורי, לא עבור הקולגות שלי. בשביל זה החלטתי להכין חבילת R פנימית.

בנה ותעד חבילה פשוטה מאוד המכילה רק כמה פונקציות גישה לנתונים המאורגנות סביב פונקציה get_snp. הכנתי גם אתר לעמיתים שלי pkgdown, כך שהם יכולים לראות בקלות דוגמאות ותיעוד.

מטמון חכם

מה למדתי: אם הנתונים שלך מוכנים היטב, האחסון במטמון יהיה קל!

מכיוון שאחד מזרימות העבודה העיקריות יישם את אותו מודל ניתוח על חבילת SNP, החלטתי להשתמש ב-binning לטובתי. בעת העברת נתונים באמצעות SNP, כל המידע מהקבוצה (bin) מצורף לאובייקט המוחזר. כלומר, שאילתות ישנות יכולות (בתיאוריה) לזרז את העיבוד של שאילתות חדשות.

# Part of get_snp()
...
  # Test if our current snp data has the desired snp.
  already_have_snp <- desired_snp %in% prev_snp_results$snps_in_bin

  if(!already_have_snp){
    # Grab info on the bin of the desired snp
    snp_results <- get_snp_bin(desired_snp)

    # Download the snp's bin data
    snp_results$bin_data <- aws.s3::s3readRDS(object = snp_results$data_loc)
  } else {
    # The previous snp data contained the right bin so just use it
    snp_results <- prev_snp_results
  }
...

בעת בניית החבילה, הרצתי אמות מידה רבות כדי להשוות מהירות בעת שימוש בשיטות שונות. אני ממליץ לא להזניח את זה, כי לפעמים התוצאות הן בלתי צפויות. לדוגמה, dplyr::filter היה הרבה יותר מהיר מאשר לכידת שורות באמצעות סינון מבוסס אינדקס, ואחזור עמודה בודדת ממסגרת נתונים מסוננת הייתה מהירה בהרבה משימוש בתחביר אינדקס.

שימו לב שהאובייקט prev_snp_results מכיל את המפתח snps_in_bin. זהו מערך של כל ה-SNPs הייחודיים בקבוצה (bin), המאפשר לך לבדוק במהירות אם כבר יש לך נתונים משאילתה קודמת. זה גם מקל על לולאה דרך כל ה-SNPs בקבוצה (פח) עם הקוד הזה:

# Get bin-mates
snps_in_bin <- my_snp_results$snps_in_bin

for(current_snp in snps_in_bin){
  my_snp_results <- get_snp(current_snp, my_snp_results)
  # Do something with results 
}

ממצאים

כעת אנו יכולים (והתחלנו להפעיל ברצינות) מודלים ותרחישים שבעבר לא היו נגישים עבורנו. הדבר הטוב ביותר הוא שעמיתיי למעבדה לא צריכים לחשוב על שום סיבוכים. יש להם פשוט פונקציה שעובדת.

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

המהירות עלתה בצורה ניכרת. אנחנו בדרך כלל סורקים שברי גנום משמעותיים מבחינה תפקודית. בעבר, לא יכולנו לעשות זאת (התברר שזה יקר מדי), אבל עכשיו, הודות למבנה הקבוצה (פח) והאחסון במטמון, בקשה ל-SNP אחד נמשכת בממוצע פחות מ-0,1 שניות, והשימוש בנתונים הוא כל כך נמוך שהעלויות עבור S3 הן בוטנים.

לאחרונה הוכנסתי לשינוי של 25+ TB של נתוני גנוטיפ גולמיים עבור המעבדה שלי. כשהתחלתי, השימוש ב-Spark לקח 8 דקות ועלה 20 דולר כדי לבצע שאילתות SNP. לאחר שימוש ב-AWK+ #rstats לעיבוד, זה לוקח כעת פחות מעשירית השנייה ועולה $10. האישיות שלי # BigData לנצח. pic.twitter.com/ANOXVGrmkk

— ניק סטרייר (@NicholasStrayer) מאי 30, 2019

מסקנה

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

כלי ביג דאטה הם מגוונים. אם יש לך זמן, כמעט בוודאות תוכל לכתוב פתרון מהיר יותר באמצעות טכניקות חכמות של ניקוי, אחסון וחילוץ נתונים. בסופו של דבר זה מסתכם בניתוח עלות-תועלת.

מה למדתי:

• אין דרך זולה לנתח 25 TB בכל פעם;
• היזהר עם גודל קבצי הפרקט והארגון שלהם;
• מחיצות ב-Spark חייבות להיות מאוזנות;
• באופן כללי, לעולם אל תנסה ליצור 2,5 מיליון מחיצות;
• המיון עדיין קשה, וכך גם הקמת Spark;
• לפעמים נתונים מיוחדים דורשים פתרונות מיוחדים;
• צבירת הניצוץ היא מהירה, אבל החלוקה עדיין יקרה;
• אל תישן כשהם מלמדים אותך את היסודות, מישהו כנראה פתר לך את הבעיה עוד בשנות השמונים;
• gnu parallel - זה דבר קסום, כולם צריכים להשתמש בו;
• Spark אוהב נתונים לא דחוסים ואינו אוהב שילוב מחיצות;
• ל-Spark יש יותר מדי תקורה בעת פתרון בעיות פשוטות;
• המערכים האסוציאטיביים של AWK יעילים מאוד;
• אתה יכול ליצור קשר stdin и stdout מסקריפט R, ולכן משתמשים בו בצנרת;
• הודות ליישום נתיב חכם, S3 יכול לעבד קבצים רבים;
• הסיבה העיקרית לבזבוז זמן היא אופטימיזציה של שיטת האחסון שלך בטרם עת;
• אל תנסה לייעל משימות באופן ידני, תן ​​למחשב לעשות זאת;
• ה-API צריך להיות פשוט למען הקלות והגמישות בשימוש;
• אם הנתונים שלך מוכנים היטב, האחסון במטמון יהיה קל!

מקור: www.habr.com