2017-07-17 20 views
1

我在FASTQ格式,DNA序列數據這需要每記錄一個4行的格式:
@序列報頭信息
序列
+
質量分數計數和除去字符在不同的行

序列行中的每個字符在質量得分線中都有相應的字符。所有序列包含8-9個鹼基的條形碼序列,接着是表示生物序列起點的引物區域。在條形碼前面,可能有或可能不存在一些垃圾箱。我需要根據條形碼將序列(及其相關信息)分成單獨的文件,並在引物之前刪除所有內容。對於序列本身,我可以使用grep做到這一點和sed(通過遍歷條形碼序列與樣本數的文件中設置a和b是變量,正則表達式打底漆):

a=AAAGCG 
b=1 
cat DNA.fastq | grep -A2 -B1 -E "${a}""CCTACGGG[ACGT]{1}GGC[AT]{1}GCAG" | grep -v "^--$" | sed -E 's/.+(CCTACGGG[ACGT]{1}GGC[AT]{1}GCAG)/\1/' > sample_${b}.fastq 

這將很好地工作在這樣的輸入數據,其中,我想從線2和ATAAAGCG的從第6行)的開始的開始除去AAAGCG:

@M00816:90:000000000-AE7TD:1:2116:19022:11483 1:N:0:1 
AAAGCGCCTACGGGTGGCAGCAGTGGGGAATTTTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGTGGGAAGAAGGCCTTCGGGTTGTAAACCACTTTTGTCAGGGAAGAAACGGTCTGAACTAATATTTCGGACTAATGACGGTACCTGAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCAGCTTTGCAAGACAGATGTGAAATCCCCGGGCTCAACCTGGGCACTGCA 
+ 
CCCCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFCFGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGFGGGGGGCFFGGGGGGGGGGGGGGGDGGGGGGGGGGGGGGGGGGGGGGGGGGCAGGGG?FGGGGGGGFGGGG7FG7CGGGGGGGGGGGGEGGGFFFG585C98DEEG?CFG*<<@>37;68CEFCFGF99EGGEEDGE54<9C<8CC>*854C<. 
@M00816:90:000000000-AE7TD:1:2118:10209:10682 1:N:0:1 
ATAAAGCGCCTACGGGGGGCAGCAGTGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTACCCGGGATGATAATGACAGTACCGGGAGAATAAGCCCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGGGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTATTCAAGTCAGAGGTGCAAGCCTGGAGCTCAACTCCAGAACTGCCTTTGAACCTGGATAGCTTGAATCAT 
+ 
[email protected]FGGGGGGEFGFGFFCFGGG,@FFFG<,<FFGGGGGGDGGGGGGGGFGG7BEE>EDGGGFGFGGFFGGGGGGGGGGGCFGGGC:CGFGGG8,<[email protected]:>G>@FGDCG=F?*1+?E>[email protected]@FFG9<[email protected]:EECC8/;;EGC?EC>FF:DA74) 

的問題是如何然後從開始去除的字符的權數的質量得分。在上面的例子中,我需要從第4行(CCCCCG)開始刪除6個字符,第8行(--CCCCCG)刪除8個字符。這可以通過計算第2行& 4中的條形碼&引物之前的鹼基數並從該分數中移除該數目的字符來完成。或者,可以在輸出文件上使用第二個命令來計算每個序列的長度(sample_1.fastq中的第2行),並將質量分數修剪爲相同的長度。兩者的問題在於它們需要存儲來自一行的信息並將其用於後續行中的過程。

理念1:存儲序列長度在單獨的文本文件

cat sample_1.fastq | sed -n '2~4p' | awk 'BEGIN{print length($0)} > seq_lengths.txt' 

的問題,這是我無法想象如何循環同時在sample_1.fastq每條記錄的文件的每一行(如與嵌套循環相反)。

理念2:

cat sample_1.fastq | awk '$0 == "+" { 
a=length $NR-1 
b=gensub(/^.{a}/, 1, $NR+1) 
print $(NR-2), $(NR-1), $0, b 
}' > sample1test.fastq 

測試這表明,長度$ NR-1給出了長度-1,而不是以前的線的長度! (我仍然是一個awk新手。)

想法3:計算sample_1.fastq中每個序列行的長度並將其打印在下一行。再次捕獲該文件,grep -B2 -A2這個數字並作爲變量存儲。使用sed將質量分數修剪爲變量的值並刪除多餘的行。這需要存儲一個變量而不會中斷通過管道的流量,根據我的理解,這是不可能的。

想法4:將2的倍數設置爲var1,將4的倍數設置爲var2。使用awk打印sample_1.fastq行號var1的長度並將其設置爲var3。將線條var2修剪爲長度爲var3。這需要能夠通過號碼呼叫線路,這我也不可能瞭解。

任何想法或見解將不勝感激!非常感謝。

+0

通過在重定向之前將文件名放在awk語句的末尾,可以避免需要使用cat進行管道操作。還需要使用長度作爲函數,即a =長度($ NR-1) –

+0

目前尚不清楚你想要做什麼。試着給出一個例子,說明你需要計算什麼,以及你必須從行中刪除什麼。 – weibeld

+0

這種格式足夠先進,我強烈建議使用更好的語言,如Python,Ruby或Java來完成這項工作。 – l0b0

回答

0

這似乎是一個不常見的需求,所以我要做的第一件事就是尋找一個預先寫好的工具。

幾個頂級的谷歌點擊:

BTW Biostars對於這樣的問題可能是一個更好的地方。