1. Bağdaştırıcı Kırpma
Bir dizi okuma ile karşılaştıktan sonra yaptığım ilk şeylerden biri, bağdaştırıcı dizilerini okumaların başından ve sonundan kaldırmaktır. Çoğu temel arama yazılımı, bir miktar yerleşik adaptör kırpma içerir, ancak neredeyse her zaman bazı adaptör sıralarının kalacağı durumdur. Bağdaştırıcıların kaldırılması eşleştirme için yararlıdır çünkü eşleştiricinin bir eşleşme yapmak için harcaması gereken çabayı azaltır (bu, bazı sınır okumalarının doğru şekilde eşleştirilmesine yardımcı olabilir).
Tercih ettiğim adaptör düzeltici Trimmomatic. Palindromik adaptörü arama yeteneğine sahiptir ve uyarlanabilir bir sürgülü pencere seçeneği içerir. Trimmomatic iş parçacıklıdır ve dosyalar arasında her seferinde bir zincir modunda çalıştırıldığında oldukça hızlı çalışıyor gibi görünüyor. İşte bazı örneklerden okumaları kırpmak için çalıştırabileceğim bir örnek komut dosyası. Bu, kırpılmış FASTQ dosyalarını bir kırpılmış alt dizine koyar:
TRIMPATH = / data / all / program / trimmomatic / Trimmomatic-0.36; JARPATH = "$ {TRIMPATH} / trimmomatic-0.36.jar "; ADAPTfile =" $ {TRIMPATH} /adapters/TruSeq3-PE-2.fa "; CMDS =" ILLUMINACLIP: $ {ADAPTfile}: 2: 30: 10: 5: true LEADING: 3 TRAILING: 3 SLIDINGWINDOW: 10: 20 MINLEN: 40 "; mkdir -p kırpılmış; INFILE1 için * _1.fastq.gz; do base = $ (basename "$ INFILE1" _1.fastq.gz); echo $ {base}; INFILE2 = "$ {temel} _2.fastq.gz"; OUTFILE_P1 = "kırpılmış / $ {baz} _P1.fastq.gz"; OUTFILE_P2 = "kırpılmış / $ {baz} _P2.fastq.gz"; OUTFILE_U1 = "kırpılmış / $ {temel} _U1.fastq.gz"; OUTFILE_U2 = "kırpılmış / $ {baz} _U2.fastq.gz"; java -jar "$ {JARPATH}" PE -iş parçacığı 20 -phred33 \ "$ {INFILE1}" "$ {INFILE2}" \ "$ {OUTFILE_P1}" "$ {OUTFILE_U1}" "$ {OUTFILE_P2}" "$ { OUTFILE_U2} "\" $ {CMDS} "; tamamlandı;
Trimmomatic kötü giriş verilerine bir şekilde duyarlıdır ve farklı uzunlukta kaliteli dizeler görürse yüksek sesle başarısız olur (yani çalışmayı durdurur) dizeleri sıralamak için:
İstisna dizisindeki "Konu-1" java.lang.RuntimeException: Dizi ve kalite uzunluğu eşleşmiyor: 'TACATGGCCCTGAAATGACTTTCACCCAGGCAACCAGTGCCCCCTGTATAGACACATGCCTTGGGCGCTCCCCACCCTTCCTCGCGTGGCCACACCTCTGT' -AAFFJFJJ7A-FJJJFJFJFJJFFA-FFAF-<A-<FFFJA-<-A7-F<<FFJJAJJJJJJJJJ' vs 7-7AJ7<AAJA - - org.usadellab.trimmomatic.fastq.FastqParser.parseOne de org.usadellab.trimmomatic.fastq.FastqRecord.<init> (FastqRecord.java:25) de J7<ACTGCTGTGGGGCACCCAGCCCCCCAGATAGCCTGGCAGAAGGATGGGGGCACAGACTTCCCAGCTGCACGGGAGAGAC'(FastqParser.java:89) --< org.usadellab.trimmomatic.fastq.FastqParser.next (FastqParser.java:179) at org.usadellab.trimmomatic.threading.ParserWorker.run (ParserWorker.java:42) java.lang.Thread.run (Thread.java) adresinde : 745)
2. Haritalamayı Okuyun
Okuma haritalama, bir hedef genom içindeki okuma için en olası konumu bulur. Son derece hızlı yaklaşık haritalayıcılar mevcuttur, ancak bunlar henüz varyant çağırma için çalışmaz ve tam bir haritalama yaklaşımı gereklidir . Şu anda tercih ettiğim eşleştiricim HISAT2. Bunu, çift okuma sorunu ve yerel varyant farkındalık eşleme nedeniyle BWA yerine kullanıyorum. HISAT2, Bowtie2 ve Tophat2 (JHUCCB) ile aynı bilgi işlem grubu tarafından yapılmıştır ve bu diğer programları değiştirmek için önerilen araçtır.
HISAT2, hem HISAT hem de TopHat2'nin halefidir. HISAT ve TopHat2 kullanıcılarının HISAT2'ye geçmelerini öneririz.
HISAT2 sayfası, SNP varyantları dahil olmak üzere insan genomu için bir genomik dizin dosyasına bir bağlantı içerir. İnsan okumalarını haritalandırırken kullanırım. HISAT2 de iş parçacıklıdır, bu yüzden onu dosyalar üzerinde teker teker çalıştırıyorum ve sıralı bir BAM dosyası oluşturmak için samtools içinden geçiyorum:
mkdir mapped for r1 in trimmed / * _ P1.fastq.gz; do base = $ (temel adı "$ {x}" _P1.fastq.gz); r2 = "kırpılmış / $ {taban} _P2.fastq.gz"; echo $ {base};
hisat2 -p 20 -t -x / data / all / genomes / hsap / hisat_grch38_snp / genome_snp -1 \ "$ {r1}" -2 "$ {r2}" 2> "eşlenmiş / hisat2 _ $ {y} _vs_grch38_stderr.txt" | \ samtools sırala > "eşlenmiş / hisat2 _ $ {y} _vs_grch38.bam"; bitti
3. Varyant Çağrısı
Bunu yapmak için şu anda samtools / bcftools kullanıyorum, ancak yakın zamanda samtools'un bir dinozor programı olduğuna ve daha iyi yaklaşımların mevcut olduğuna dair bir hibe gözden geçiren yanıtı aldığımız için diğer fikirlerle ilgilenebilirim. Samtools şu anda varyant çağrısı için evreler ile çalışmıyor, bu yüzden komutları bir metin dosyasına kaydedip GNU Parallel üzerinden çalıştırıyorum. Bu, genom dosyalarının ilk önce samtools faidx aracılığıyla indirilmesini ve dizine eklenmesini gerektirir. Bu, varyantlar dizininde bir dizi gzipli VCF dosyası oluşturur:
mkdir -p varyantları; (eşlenmiş / *. Bam içindeki x için; do echo samtools mpileup - v -f /data/all/genomes/hsap/hisat_grch38_snp/Homo_sapiens.GRCh38.dna.primary_assembly.fa "$ {x}" \ | \ bcftools call --ploidy GRCh38 -v -m -O z -o varyantları / $ (temel adı "$ {x}" .bam) .vcf.gz; tamamlandı) > call_jobs.txtcat call_jobs.txt | paralel -u;
Yavaş, ancak daha doğru bir alternatif olarak, tüm örnekler için aynı anda varyantlar çağrılabilir:
samtools mpileup -v -f /data/all/genomes/hsap/hisat_grch38_snp/Homo_sapiens.GRCh38.dna.primary_assembly.fa eşlenmiş / *. bam | \ bcftools call --ploidy GRCh38 -v -m -O z -o varyantları / hisat2_allCalled_vs_grch38.vcf.gz
Daha hızlı bir işlem isteniyorsa, varyant çağırma kromozomlar arasında paralelleştirilebilir paralel kullanarak ve daha sonra bcftools norm kullanılarak birleştirildi. Bunun uygulanması okuyucuya bir alıştırma olarak bırakılmıştır.