使用python讀取和分析fasta檔案

在生物資訊學中，fasta格式（又稱為pearson格式）是一種基於文字的、用於表示核苷酸序列或氨基酸序列的格式。fasta檔案以序列標識和序列作為乙個基本單元，每個基本單元分為兩部分：序列標記和序列本身。第一行以『>』開頭，後面緊跟序列標記；從第二行開始，直到下乙個標記行(『>』開頭行)出現，或檔案末尾，這部分為序列本身。值得注意的是，序列中的換行符應該被忽略。

fasta檔案格式詳細介紹請參見：傳送門1

傳送門2

例如：

>seq1
atcgatcg
>seq2
aaatttcc
cggg
>seq3
tagctagctagc

上面就是乙個fasta檔案的基本格式，描述了3條dna序列seq1、seq2和seq3。值得注意的是fasta檔案允許在序列中存在換行，這給我們的分析工作造成了麻煩。

perl讀取fasta檔案

編寫perl指令碼讀取fasta比較簡單，我們可以更改perl的特殊變數$/，讓"行"的概念由換行符分割改為由fasta檔案中序列標記開頭標識">"分割。請參考以下**

open in,"<",$fa or die $!;
#$/ 預設為換行符'\n'，我們將其改為'>'
$/ = ">";;
while()
#為了避免對之後**的影響，及時把$/改回來
$/ = "\n";
close in;

我封裝了一些處理fasta檔案的函式，涉及到讀取序列、滑窗操作等，在這裡分享給大家。

def readfa(fa):
'''@msg: 讀取乙個fasta檔案
@param fa   fasta 檔案路徑
@return:  返回乙個生成器，能迭代得到fasta檔案的每乙個序列名和序列
'''with open(fa,'r') as fa:
seqname,seq='',''
while 1:
line=fa.readline()
line=line.strip('\n')
if (line.startswith('>') or not line) and seqname:
yield((seqname,seq))
if line.startswith('>'):
seqname = line[1:]
seq=''
else:
seq+=line
if not line:break
def getseq(fa,queryseqname,start=1,end=0):
'''@msg: 獲取fasta檔案的某一條序列
@param fa   fasta 檔案路徑
@param queryseqname   序列名
@param start   擷取該序列時，起始位置，可省略，預設為1
@param end   fasta 擷取該序列時，最後位置，可省略，預設為該序列全長
@return:  返回找到(擷取到)的序列
'''if start<0: start=start+1
for seqname,seq in readfa(fa):
if queryseqname==seqname:
if end!=0: returnseq = seq[start-1:end];print(start-1)
else: returnseq = seq[start-1:]
return returnseq
def getreversecomplement(sequence):
'''@msg: 獲取反向互補序列
@param sequence   一段dna序列
@return:  返回反向互補序列
'''sequence = sequence.upper()
sequence = sequence.replace('a', 't')
sequence = sequence.replace('t', 'a')
sequence = sequence.replace('c', 'g')
sequence = sequence.replace('g', 'c')
return sequence.upper()[::-1]
def getgc(sequence):
'''@msg: 獲取某一條序列的gc含量
@param sequence   一段dna序列
@return:  返回gc含量
'''sequence=sequence.upper()
return (sequence.count("g")+sequence.count("c"))/len(sequence)
def readseqbywindow(sequence,winsize,stepsize):
'''@msg: 滑窗讀取某一條序列
@param sequence   一段dna序列
@param winsize   視窗大小
@param stepsize   步長
@return:   返回乙個生成器，可迭代得到該序列的每乙個視窗序列
'''if stepsize<=0: return false
now = 0
seqlen = len(sequence)
while(now+winsize-stepsize以下是兩個小例子：
fa="./sequence.fasta"
#讀取乙個fasta檔案，並輸出其中的每一條序列名，序列長度和gc含量
for seqname,seq in readfa(fa):
seqlen = len(seq)
gc = getgc(seq)
print(seqname,seqlen,gc)
#讀取乙個fasta檔案，並以1000bp為視窗、100bp為步長讀取每條序列，
#計算每個視窗的gc含量，並記錄在字典中
gclst = {}
for seqname,seq in readfa(fa):
gclst.setdefault(seqname,)
for winseq in readseqbywindow(seq,1000,100):

使用python讀取和分析fasta檔案

Python使用pandas讀取csv和excel

python讀取使用json

使用Python讀取幾何

使用python讀取和分析fasta檔案

Python使用pandas讀取csv和excel

python讀取使用json

使用Python讀取幾何

相關推薦