搜索
查看: 10572|回复: 41

生信编程直播第0题-生信编程很简单!

[复制链接]

633

主题

1172

帖子

3947

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3947
发表于 2017-2-8 08:09:23 | 显示全部楼层 |阅读模式
有不少学员反映前面的题目有点太难了,我们临时插播这个系列!
分成3次录制!

有很多人认为生信工程师可以不需要学习编程,这一点,我反对!
有很多人认为学习编程来解决生物信息数据处理很难,这一点,我也反对!
首先安装bowtie2这个软件,里面有测试数据!
  • ## Download and install bowtie
  • cd ~/biosoft
  • mkdir bowtie &&  cd bowtie
  • wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.2.9/bowtie2-2.2.9-linux-x86_64.zip  
  • unzip bowtie2-2.2.9-linux-x86_64.zip

针对测试数据,我会根据下面列出的需求现成编程讲解!
对FASTQ的操作
  • 5,3段截掉几个碱基
  • 序列长度分布统计
  • FASTQ 转换成 FASTA
  • 统计碱基个数及GC%
对FASTA的操作
  • 取互补序列
  • 取反向序列
  • DNA to RNA
  • 大小写字母形式输出
  • 每行指定长度输出序列
  • 按照序列长度/名字排序
  • 提取指定ID的序列
  • 随机抽取序列
高级难度:






上一篇:生信人的windows电脑必须安装的软件
下一篇:对数据做对数转换会改变相关性吗?
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

11

主题

45

帖子

232

积分

中级会员

Rank: 3Rank: 3

积分
232
发表于 2017-3-18 19:04:01 | 显示全部楼层
生信编程直播第零题第一讲,Fastq文件简单统计处理
by dongye

浓缩的精华-让你的代码高大上 视频代码
[Python] 纯文本查看 复制代码
# /bin/python
# encoding=utf-8


"""
简介: fastq文件数据处理工具
作者: dongye
时间: 2017年3月18日01:11:16
"""


from optparse import OptionParser
import sys

args = sys.argv


def cut_fastq(filename, cut_5=0, cut_3=0):
    """
    从5'端,3'端截掉几个碱基

    :param filename: fastq文件名
    :param cut_5: 5'端截掉碱基数, 大于0
    :param cut_3: 3'端截掉碱基数, 大于0
    :return:
    """

    if cut_3 < 0:
        raise ValueError('cut_3: %s < 0' % cut_3)
    if cut_5 < 0:
        raise ValueError('cut_5: %s < 0' % cut_5)

    with open(filename) as if_fq:
        while if_fq:
            line_1 = if_fq.readline().strip('\n')
            if not (if_fq and line_1):
                break
            line_2 = if_fq.readline().strip('\n')
            line_3 = if_fq.readline().strip('\n')
            line_4 = if_fq.readline().strip('\n')

            print line_1
            print cut_3 and line_2[cut_5:][:-cut_3] or line_2[cut_5:]
            print line_3
            print cut_3 and line_4[cut_5:][:-cut_3] or line_4[cut_5:]


def length_distribution(filename):
    """
    统计fastq文件reads长度分布

    :param filename: fastq文件名
    :return:
    """

    reads_len = {}
    with open(filename) as if_fq:
        while if_fq:
            line_1 = if_fq.readline().strip('\n')
            if not (if_fq and line_1):
                break
            line_2 = if_fq.readline().strip('\n')
            line_3 = if_fq.readline().strip('\n')
            line_4 = if_fq.readline().strip('\n')

            length = len(line_2)
            if not length in reads_len:
                reads_len[length] = 0
            reads_len[length] += 1

    for length, num in reads_len.items():
        print "%s\t%s" % (length, num)
    return reads_len


def fq_2_fa(filename):
    """
    将fastq格式文件转换为fasta格式

    :param filename: fastq文件
    :return:
    """

    with open(filename) as if_fq:
        while if_fq:
            line_1 = if_fq.readline().strip('\n')
            if not (if_fq and line_1):
                break
            line_2 = if_fq.readline().strip('\n')
            line_3 = if_fq.readline().strip('\n')
            line_4 = if_fq.readline().strip('\n')

            print '>'+line_1
            print line_2


def fq_length(filename):
    """
    统计fastq文件中碱基数量(reads总长度)

    :param filename: fastq文件
    :return:
    """

    num = 0
    with open(filename) as if_fq:
        while if_fq:
            line_1 = if_fq.readline().strip('\n')
            if not (if_fq and line_1):
                break
            line_2 = if_fq.readline().strip('\n')
            line_3 = if_fq.readline().strip('\n')
            line_4 = if_fq.readline().strip('\n')

            num += len(line_2)
    print num
    return num


def count_gc(filename):
    """
    统计fastq文件中reads的GC含量

    :param filename: fastq文件
    :return:
    """

    GC = 0
    with open(filename) as if_fq:
        while if_fq:
            line_1 = if_fq.readline().strip('\n')
            if not (if_fq and line_1):
                break
            line_2 = if_fq.readline().strip('\n').upper()
            line_3 = if_fq.readline().strip('\n')
            line_4 = if_fq.readline().strip('\n')
            GC += line_2.count('G') + line_2.count('C')

    # 函数的好处
    SUM = fq_length(filename)
    print "GC : %s" % (GC * 1.0 / SUM)
    return GC * 1.0 / SUM


def main(args):
    parser = OptionParser()
    parser.add_option("-f", "--fastq", dest="filename",
                      help="fastq filename", metavar="FILE")

    parser.add_option("--cut", "--cut-fastq", dest="cut",
                      help="cut fastq file", action="store_true",
                      default=False)
    parser.add_option("-5", "--cut-5", dest="cut_5", type=int,
                      help="cut fastq N(N>0) bases from 5'", metavar="INT", default=0)
    parser.add_option("-3", "--cut-3", dest="cut_3", type=int,
                      help="cut fastq N(N>0) bases from 3'", metavar="INT", default=0)

    parser.add_option("--len-dis", dest="len_dis",
                      help="sequence length distribution", action="store_true",
                      default=False)

    parser.add_option("--fasta", "--fastq2fasta", dest="fasta",
                      help="convert fastq to fasta", action="store_true",
                      default=False)

    parser.add_option("--len", "--fastq-length", dest="fq_len",
                      help="total reads number", action="store_true",
                      default=False)

    parser.add_option("--gc", "--count-gc", dest="gc",
                      help="print GC%", action="store_true",
                      default=False)

    (options, args) = parser.parse_args()
    if not options.filename:
        parser.print_help()
    filename = options.filename

    if options.cut:
        cut_5 = options.cut_5
        cut_3 = options.cut_3
        cut_fastq(filename, cut_5, cut_3)

    if options.len_dis:
        length_distribution(filename)

    if options.fq_len:
        fq_length(filename)

    if options.fasta:
        fq_2_fa(filename)

    if options.gc:
        count_gc(filename)


if __name__ == '__main__':
    main(args)


回复 支持 3 反对 0

使用道具 举报

0

主题

5

帖子

89

积分

注册会员

Rank: 2

积分
89
发表于 2017-4-8 20:12:07 | 显示全部楼层
本帖最后由 lingboling002 于 2017-4-8 20:17 编辑

生信第二次作业,fastq文件的简单处理。大学自学自考的c语言二级三级,后来就扔掉了,这次正好重新复习。随机拿了数据进行练习,基本信息如图。
序列长度分布统计
                                                      
[Perl] 纯文本查看 复制代码
perl -alne '{print}' 39.fastq  |head #先看fastq文件结构(每4行一个reads)
perl -alne '{print if $.%4==1}' 39.fastq  |head #看每个reads的第一行
perl -alne '{print length }'  39.fastq |head #看每行长度
perl -alne '{print length if $.%4==2  }'  39.fastq |head #再看每个reads的序列行长度(数据质量较好,大部分长度在150)
perl -alne '{print length if $.%4==2  }'  39.fastq >length.csv #将上一步序列长度数据输出为csv文件(用于后续验证)
perl -alne '{print length if $.%4==2  }'  39.fastq |perl -alne '{$h{$_}++  }END{print "$_\t$h{$_}" foreach sort{$a <=> $b} keys %h}'  |head #求频数
grep 35 length.csv |wc #查看csv中长度包含35的总和(伪)(若已知序列长度最大值为三位数,可以再计算长度为135的总和,再求差值

结果:
39.fastq文件共有14040924行,为reads的4倍。
hash输出长度为35的序列有109个。csv文件中grep 35之后有5602个结果,grep 135之后有5493个结果,差值即为109。
通过R将csv文件进行作图如下

fastq文件转fasta文件

删除fastq中每个reads的第3、4行,再将第1行的@改为>。
[Perl] 纯文本查看 复制代码
perl -alne '{print if ($.%4==1 | $.%4==2)}' 39.fastq | perl -alne '{s/@/>/g; print}' >fasta.csv

结果
>E00491:27:H37GGALXX:6:1101:17401:1924 1:N:0:GGCTAC
TTCTATAGGAACTGCCAAGAAAACAGGCGATCAGGAGAGCAGAGGAAGTGGGGATATAGGTGCTATCAGTGGGGGCAGGGAGCTATGCGCTGGCAGGGCCGCCGCTGAGCGTGTGAAACACAAGTCCTAGAGGAGGAAGAAACTACCACT
>E00491:27:H37GGALXX:6:1101:14884:1959 1:N:0:GGCTAC
CGGTTTTTGTATATAAAATTCATGTTTCCAATCTCTCTCTCCTTGGTCGGTGACACTGACTAGCCTATCTTGAACAGATATTTATTTTTTCTAACACTCAGCTCTGCGCTCCCCGATTCCCTGTCTCGTCTGCACGATTTGCATTGAAAT







本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 1 反对 0

使用道具 举报

2

主题

34

帖子

707

积分

高级会员

Rank: 4

积分
707
发表于 2017-3-27 19:45:54 | 显示全部楼层
本帖最后由 x2yline 于 2017-4-3 09:00 编辑

077-x2yline
(1)处理fastq文件python3代码如下:
[Python] 纯文本查看 复制代码
# encoding=utf-8
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
import os
 
def density_plot(data_list,xlabel,title):
    if max(data_list) == min(data_list):
        print('\nAll data is same, can not plot the density line!!!!!!\n')
    elif len(data_list) > 500000:
        x0 = []
        y0 = []
        for i in sorted(list(set(data_list))):
            x0.append(i)
            y0.append(data_list.count(i))
        fig = plt.figure(1)
        fig.patch.set_facecolor('w')
        ax = fig.add_axes([0.2,0.2,0.5,0.6])
        data_range = max(data_list) - min(data_list)
        ax.plot(x0, y0,'b*')
        # axis seting
        ax.set_xlim([int(min(data_list)-data_range/10),int(max(data_list)+data_range/10)])
        ax.spines['top'].set_visible(False)
        ax.spines['right'].set_visible(False)
        ax.get_xaxis().tick_bottom()
        ax.get_yaxis().tick_left()
        ax.get_xaxis().set_tick_params(direction='out')
        ax.get_yaxis().set_tick_params(direction='out')
        ax.spines['left'].set_position(('outward', 10))
        ax.spines['bottom'].set_position(('outward', 10))
        ax.set_xlabel(xlabel)
        ax.set_xticks(np.linspace(int(min(data_list)-data_range/10),int(max(data_list)+data_range/10),6))
        # text seting
        ax_stats = fig.add_axes([0.76,0.2,0.2,0.6])
        ax_stats.set_axis_off()
        ax_stats.text(0, 0.8,'Mean: {:.3f}'.format(float(np.mean(data_list))), fontsize=12, color='blue')
        ax_stats.text(0, 0.6,'Max: {:.3f}'.format(max(data_list)), fontsize=12, color='red')
        ax_stats.text(0, 0.4,'Min: {:.3f}'.format(min(data_list)), fontsize=12, color='green')
        ax_stats.text(0, 0.2,'Median: {:.3f}'.format(float(np.median(data_list))), fontsize=12,color='blue')
        fig.suptitle('\n\n\n'+title)
        fig.savefig('density_'+xlabel, dpi=150)
        plt.show()
        print('\nFigure saved in {0}\n'.format(str(os.path.join(os.getcwd(),'density_'+xlabel))))
    else:
        fig = plt.figure(1)
        fig.patch.set_facecolor('w')
        ax = fig.add_axes([0.2,0.2,0.5,0.6])
        density = stats.kde.gaussian_kde(data_list)
        data_range = max(data_list) - min(data_list)
        x = np.arange(min(data_list),max(data_list),data_range/40)
        ax.plot(x, density(x),'b-')
        # axis seting
        ax.set_xlim([int(min(data_list)-data_range/10),int(max(data_list)+data_range/10)])
        ax.spines['top'].set_visible(False)
        ax.spines['right'].set_visible(False)
        ax.get_xaxis().tick_bottom()
        ax.get_yaxis().tick_left()
        ax.get_xaxis().set_tick_params(direction='out')
        ax.get_yaxis().set_tick_params(direction='out')
        ax.spines['left'].set_position(('outward', 10))
        ax.spines['bottom'].set_position(('outward', 10))
        ax.set_xlabel(xlabel)
        ax.set_xticks(np.linspace(int(min(data_list)-data_range/10),int(max(data_list)+data_range/10),6))
        # text seting
        ax_stats = fig.add_axes([0.76,0.2,0.2,0.6])
        ax_stats.set_axis_off()
        ax_stats.text(0, 0.8,'Mean: {:.3f}'.format(float(np.mean(data_list))), fontsize=12, color='blue')
        ax_stats.text(0, 0.6,'Max: {:.3f}'.format(max(data_list)), fontsize=12, color='red')
        ax_stats.text(0, 0.4,'Min: {:.3f}'.format(min(data_list)), fontsize=12, color='green')
        ax_stats.text(0, 0.2,'Median: {:.3f}'.format(float(np.median(data_list))), fontsize=12,color='blue')
        fig.suptitle('\n\n\n'+title)
        fig.savefig('density_'+xlabel, dpi=150)
        plt.show()
        print('\nFigure saved in {0}\n'.format(str(os.path.join(os.getcwd(),'density_'+xlabel))))
     
def cut_53(file, cut_5, cut_3,buffer=2048*3096):
    with open(file) as f:
    	name = f.read(4).strip().split('@')[-1]
    with open(file, 'r') as f:
        head = ''
        while True:
            raw_tmp = f.read(buffer)
            tmp = (head + raw_tmp).split('\n@' + name)
            tmp_records = tmp[:-1]
            head = tmp[-1]
            for i in tmp_records:
                if not i.startswith('@'):
                    i = '@' + name + i
                if i:
                    i = i.split('\n')
                    print(i[0])
                    if len(i[1])-cut_5 <= cut_3:
                        raise ValueError('cut too much')
                    print(i[1][cut_5:-(cut_3)])
                    print(i[2])
                    print(i[3][cut_5:-(cut_3)])
            if not raw_tmp:
                break
        if head:
            if not head.startswith('@'):
                i = '@' +name+ head
            else:
                i = head
            i = i.split('\n')
            print(i[0])
            if len(i[1])-cut_5 <= cut_3:
                raise ValueError('cut too much')
            print(i[1][cut_5:-(cut_3)])
            print(i[2])
            print(i[3][cut_5:-(cut_3)])
    return(0)
                 
 
def fastq2fasta(file, buffer=2048*3069):
    head = ''
    with open(file) as f:
        name = f.read(4).strip().split('@')[-1]
    with open(file,'r') as f:
        while True:
            raw_tmp = f.read(buffer)
            tmp = (head + raw_tmp).split('\n@'+name)
            tmp_records = tmp[:-1]
            head = tmp[-1]
            for i in tmp_records:
                if i:
                    if i.strip().startswith('@'):
                        i = i.strip()[1:]
                    else:
                        i = name + i.strip()
                    fasta = '>'+'\n'.join(i.strip().split('\n')[0:2])
                    print(fasta)
            if not raw_tmp:
                break
        if head:
            for i in head.split('\n@'+name):
                if i:
                    if i.strip().startswith('@'):
                        i = i.strip()[1:]
                    else:
                        i = name + i.strip()
                    fasta = '>'+'\n'.join(i.split('\n')[0:2])
                    print(fasta)
 
def length_count(file, buffer=2400*2400):
    length_list = []
    with open(file) as p:
        name = p.read(4).strip().split('@')[-1]

    with open(file) as f:
        head = ''
        base_num = 0
        while True:
            a = f.read(buffer)
            raw_list = (head + a).split('\n@'+name)
            deal_list = raw_list[:-1]
            head = raw_list[-1]
            for record in deal_list:
                if record:
                    sequence = record.split('\n')[1].upper().strip()
                    length_list.append(len(sequence))
                    base_num += len(sequence)
            if not len(a):
                break
        if head:
            record = head
            sequence = record.split('\n')[1].upper().strip()
            length_list.append(len(sequence))
            base_num += len(sequence)
    #print(len(length_list))

    max_len = max(length_list)
    min_len = min(length_list)
    if not min_len ==max_len:
        try:
            pass
            #density_plot(length_list, xlabel='length',title='Distribution of length')
        except:
            pass
    result = ('\nMean length is {0}\nmaxium length is {1}\nminium length is {2}\n'.format(sum(length_list)/len(length_list),max_len,min_len))
    print(result)
    return(result)
 
 
def GC_count(file, buffer=2400*3200):
    with open(file) as f:
        name = f.read(4).strip().split('@')[-1]
    with open(file) as f:
        head = ''
        gc_list = []
        gc_num = 0
        base_num = 0
        while True:
            a = f.read(buffer)
            raw_list = (head + a).split('\n@' + name)
            deal_list = raw_list[:-1]
            head = raw_list[-1]
            for record in deal_list:
                if record:
                    sequence = record.split('\n')[1].upper().strip()
                    gc_item = sequence.count("G")+sequence.count("C")
                    gc_list.append(gc_item/len(sequence))
                    gc_num += gc_item
                    base_num += len(sequence)
            if not len(a):
                break
        if head:
            record = head
            sequence = record.split('\n')[1].upper().strip()
            gc_item = sequence.count("G")+sequence.count("C")
            gc_list.append(gc_item/len(sequence))
            gc_num += gc_item
            base_num += len(sequence)
    print("\nThe GC count is {0}\nThe number of bases is {1}\nThe number of sequences is {2}".format(gc_num/base_num,base_num,len(gc_list)))
    #print("Ploting...\n")
    #density_plot(gc_list, xlabel="Mean GC content",title="GC distribution over all sequences")
    return(0)
 
def main():
    import argparse
    parser = argparse.ArgumentParser()
    parser.add_argument('file', help = 'The fastq file_path')
    parser.add_argument('-c3','--cut3',type=int, help="the cut bases in the 3' of the sequences")
    parser.add_argument('-c5','--cut5',type=int,help = "the cut bases in the 5' of the sequences")
    parser.add_argument('-f', '--fasta',action="store_true",help="fastq to fasta")
    parser.add_argument('-gc', '--count',action="store_true",help = "count GC and the number of bases of fastq file")
    parser.add_argument('-l','--length',action="store_true", help="count the length of fastq file")
    args = parser.parse_args()
    if args.fasta:
        fastq2fasta(args.file)
    if args.count:
        #print("Waiting for about 20 seconds...\n")
        GC_count(args.file)

    if args.length:
        length_count(args.file)
    cut3 = 0
    cut5 = 0
    if args.cut3:
        if args.cut3 >= 0:
            cut3 =args.cut3
        else:
            raise ValueError('cut value can not be negative!!!')
    if args.cut5:
        if args.cut5 >= 0:
            cut5 = args.cut5
        else:
            raise ValueError('cut value can not be negative!!!')
    if cut3 + cut5 != 0:
        cut_53(args.file, cut5,cut3)
 
if __name__ == '__main__':
    import time
    time1 = time.time()
    main()
    time2 = time.time()
    #print('Used %s s'%(time2-time1))

直接保存文fastqc_analysis.py
命令行下运行(gc和长度统计都进行了可视化处理,需安装matplotlib,numpy和scipy,我用的示例数据是一个RNA的测序数据SRR3336676_1.fastq):
python fastqc_analysis.py SRR3336676_1.fastq -gc
gc 含量统计结果:
The GC count is 0.4646087852875249
The number of bases is 154156958
The number of sequences is
2094667


或者
python fastqc_analysis.py SRR3336676_1.fastq -l

长度统计结果:
mean length is 73.59497142027826
maxium length is 75
minium length is 32


或者
python fastqc_analysis.py SRR3336676_1.fastq -c3 5 -c5 4 > new.fastq
3’端去除5个,5’端去除4个碱基
或者
python fastqc_analysis.py SRR3336676_1.fastq -f > fastq_2fasta.fa
结果就不贴出来了,太长了

fastqc软件分析示例数据的结果与python代码的结果相同:



看完老师的视频,感觉对于经常可以使用的程序,写能够重复使用的代码才是最明智的选择,先贴上代码,后面再来补充。。。

(2)处理fasta文件的python3代码如下:
[Python] 纯文本查看 复制代码
def complementary_seq(file, buffer=3096*2048):
    base_list = ['A','T','C','G','U','N']
    base_comp_list = ['T','A','G','C','A','N']
    comp_dict = dict(zip(base_list, base_comp_list))
    head = ''
    with open(file, 'r') as f:
        while True:
            raw_tmp = f.read(buffer)
            tmp = (head + raw_tmp).split('\n')
            tmp_records = tmp[:-1]
            head = tmp[-1]
            for i in tmp_records:
                if i:
                    if i.startswith('>'):
                        print(i, end = '\n')
                    else:
                        new_seq = ''
                        for j in i.upper():
                            new_seq += comp_dict[j]
                        print(new_seq, end='\n')
            if not raw_tmp:
                break
        if head:
            if head.startswith('>'):
                print(i, end = '\n')
            else:
                new_seq = ''
                for j in head.upper():
                    new_seq += comp_dict[j]
                print(new_seq, end='\n')
    return(0)

def reverse_seq(file, buffer = 3096*2048):
    head = ''
    with open(file, 'r') as f:
        while True:
            raw_tmp = f.read(buffer)
            tmp = (head + raw_tmp)
            if tmp.count('\n>') >=2:
                tmp_list = tmp.split('\n>')
                head = '\n>' + tmp_list[-1]
                tmp_record = tmp_list[:-1]
                for i in tmp_record:
                    if not i:
                        continue
                    seq_head = i.split('\n')[0]
                    seq = ''.join(i.split('\n')[1:])
                    print(seq_head)
                    print(seq[::-1])
            else:
                head = tmp
            if not raw_tmp:
                break
        if head:
                tmp_list = head.split('\n>')
                tmp_record = tmp_list
                for i in tmp_record:
                    if not i:
                        continue
                    seq_head = i.split('\n')[0]
                    seq = ''.join(i.split('\n')[1:])
                    print(seq_head)
                    print(seq[::-1])
    return(0)

def DNA_RNA(file, buffer = 3096*2048):
    base_list = ['T','U']
    base_comp_list = ['U','T']
    comp_dict = dict(zip(base_list, base_comp_list))
    with open(file, 'r') as f:
        head = ''
        while True:
            raw_tmp = f.read(buffer)
            tmp = (head + raw_tmp).split('\n')
            tmp_records = tmp[:-1]
            head = tmp[-1]
            for i in tmp_records:
                if i:
                    if i.startswith('>'):
                        print(i, end = '\n')
                    else:
                        if 'U' in i.upper():
                            j = 'U'
                        else:
                            j = 'T'
                        new_seq = i.upper().replace(j,comp_dict[j])
                        print(new_seq, end='\n')
            if not raw_tmp:
                break
        if head:
            if head.startswith('>'):
                print(i, end = '\n')
            else:
                if 'U' in i.upper():
                    j = 'U'
                else:
                    j = 'T'
                new_seq = i.upper().replace(j,comp_dict[j])
                print(new_seq, end='\n')
    return(0)

def ToUpper(file, buffer=3069*2048):
    with open(file, 'r') as f:
        head = ''
        while True:
            raw_tmp = f.read(buffer)
            tmp = (head + raw_tmp).split('\n')
            tmp_records = tmp[:-1]
            head = tmp[-1]
            for i in tmp_records:
                if i:
                    if i.startswith('>'):
                        print(i, end = '\n')
                    else:
                        print(i.upper(), end='\n')
            if not raw_tmp:
                break
        if head:
            if head.startswith('>'):
                print(i, end = '\n')
            else:
                print(head.upper(), end='\n')
    return(0)

def ToLower(file, buffer=3069*2048):
    with open(file, 'r') as f:
        head = ''
        while True:
            raw_tmp = f.read(buffer)
            tmp = (head + raw_tmp).split('\n')
            tmp_records = tmp[:-1]
            head = tmp[-1]
            for i in tmp_records:
                if i:
                    if i.startswith('>'):
                        print(i, end = '\n')
                    else:
                        print(i.lower(), end='\n')
            if not raw_tmp:
                break
        if head:
            if head.startswith('>'):
                print(i, end = '\n')
            else:
                print(head.lower(), end='\n')
    return(0)

def OutPutLength(file, length, buffer=2048*3096):
    if length > 0:
        with open(file, 'r') as f:
            head = ''
            while True:
                raw_tmp = f.read(buffer)
                tmp = (head + raw_tmp).split('\n')
                tmp_records = tmp[:-1]
                head = tmp[-1]
                for i in tmp_records:
                    if i:
                        if i.startswith('>'):
                            print(i, end = '\n')
                        else:
                            try:
                                print(i[:length], end='\n')
                            except:
                                print(i)
                if not raw_tmp:
                    break
            if head:
                if head.startswith('>'):
                    print(i, end = '\n')
                else:
                    print(head[:length], end='\n')
    else:
        with open(file, 'r') as f:
            head = ''
            while True:
                raw_tmp = f.read(buffer)
                tmp = (head + raw_tmp).split('\n')
                tmp_records = tmp[:-1]
                head = tmp[-1]
                for i in tmp_records:
                    if i:
                        if i.startswith('>'):
                            print(i, end = '\n')
                        else:
                            try:
                                print(i[length:], end='\n')
                            except:
                                print(i)
                if not raw_tmp:
                    break
            if head:
                if head.startswith('>'):
                    print(i, end = '\n')
                else:
                    print(head[length:], end='\n')
    return(0)

def seq_order(file, TYPE, buffer=2048*3096):
    fasta_dict = {}
    with open(file, 'r') as f:
        head = ''
        while True:
            raw_tmp = f.read(buffer)
            tmp = (head + raw_tmp).split('\n')
            tmp_records = tmp[:-1]
            head = tmp[-1]
            for i in tmp_records:
                if i:
                    if i.startswith('>'):
                        KEY = i.strip()
                        fasta_dict[KEY] = ''
                    else:
                        fasta_dict[KEY] += i.strip()
            if not raw_tmp:
                break
        if head:
            if head.startswith('>'):
                KEY = head.strip()
                fasta_dict[KEY] = ''
            else:
                fasta_dict[KEY] += head.strip()
    if TYPE.lower() == 'l+':
        SORT = sorted(fasta_dict.items(), key=lambda e:len(e[1]), reverse=False)
        for item in SORT:
            print(item[0]+'\n'+item[1])
    elif TYPE.lower() == 'l-':
        SORT = sorted(fasta_dict.items(), key=lambda e:len(e[1]), reverse=True)
        for item in SORT:
            print(item[0]+'\n'+item[1])
    elif TYPE.lower() == 'n+':
        SORT = sorted(fasta_dict.items(), key=lambda e:e[0][1:], reverse=False)
        for item in SORT:
            print(item[0]+'\n'+item[1])
    elif TYPE.lower() == 'n-':
        SORT = sorted(fasta_dict.items(), key=lambda e:e[0][1:], reverse=True)
        for item in SORT:
            print(item[0]+'\n'+item[1])
    return(0)

def extract_from_id(file, ID, buffer=2048*3096):
    next_id = False
    found = False
    if ID.startswith('>'):
        ID = ID[1:]
    with open(file, 'r') as f:
        head = ''
        while True:
            raw_tmp = f.read(buffer)
            tmp = (head + raw_tmp).split('\n')
            tmp_records = tmp[:-1]
            head = tmp[-1]
            tmp_text = '\n'.join(tmp_records)
            if ID+'\n' in tmp_text:
                found = True
                if '>' in tmp_text[tmp_text.find(ID):]:
                    next_id = True
                    break
                while True:
                    raw_tmp = f.read(buffer)
                    tmp = (head + raw_tmp).split('\n')
                    tmp_records = tmp[:-1]
                    head = tmp[-1]
                    if '>' in  '\n'.join(tmp_records):
                        next_id = True
                        tmp_text += '\n'.join(tmp_records).split('>')[0]
                        break
                    else:
                        tmp_text += '\n'.join(tmp_records)
                    if not raw_tmp:
                        break
                break
            if not raw_tmp:
                break
    if found:
        all_text = tmp_text[tmp_text.find(ID):]
        if not next_id:
            all_text += '\n' + head
        if '>' in all_text:
            print('>'+all_text[:all_text.find('>')])

        else:
            print('>' + all_text)
    else:
        print('Not found!!!!!!!!!!!!!!!!!!!!!!!!!')

    return(0)

def random_seq(file, buffer=2048*3096):
    fasta_dict = {}
    with open(file, 'r') as f:
        head = ''
        while True:
            raw_tmp = f.read(buffer)
            tmp = (head + raw_tmp).split('\n')
            tmp_records = tmp[:-1]
            head = tmp[-1]
            for i in tmp_records:
                if i:
                    if i.startswith('>'):
                        KEY = i.strip()
                        fasta_dict[KEY] = ''
                    else:
                        fasta_dict[KEY] += i.strip()
            if not raw_tmp:
                break
        if head:
            if head.startswith('>'):
                KEY = head.strip()
                fasta_dict[KEY] = ''
            else:
                fasta_dict[KEY] += head.strip()
    import random
    KEY = random.choice(list(fasta_dict.keys()))
    print(KEY)
    print(fasta_dict[KEY])
    return(0)

def main():
    import argparse
    parser = argparse.ArgumentParser()
    parser.add_argument('file', help = 'The fastq file_path')
    parser.add_argument('-c', '--comp',action="store_true",help="fasta to complementary fasta")
    parser.add_argument('-r', '--reverse', action="store_true",help="fasta to reverse fasta")
    parser.add_argument('-d', '--D_R', action="store_true",help="switch between DNA and RNA fasta")
    parser.add_argument('-u', '--upper', action='store_true',help='fasta to upper fasta')
    parser.add_argument('-l', '--lower', action='store_true',help='fasta to lower fasta')
    parser.add_argument('-s', '--selection', type=int, help='select length to output, length can either be positive or negative')
    parser.add_argument('-o', '--order', help ='order the fasta file by name for n+ or n-, order the fasta file by length for l+ or l-')
    parser.add_argument('-f', '--find', help = 'find the sequence for a fasta file')
    parser.add_argument('-m', '--randm', action='store_true', help = 'extract a random sequence from a fasta file')
    args = parser.parse_args()

    if args.reverse:
        reverse_seq(args.file)
    if args.comp:
        complementary_seq(args.file)
    if args.D_R:
        DNA_RNA(args.file)
    if args.lower:
        ToLower(args.file)
    if args.upper:
        ToUpper(args.file)
    if args.selection:
        OutPutLength(args.file, args.selection, buffer=2048*3096)
    if args.order:
        seq_order(args.file, args.order, buffer=2048*3096)
    if args.find:
        extract_from_id(args.file, args.find, buffer=2048*3096)
    if args.randm:
        random_seq(args.file, buffer=2048*3096)
if __name__ == '__main__':
    main()

上述代码均使用read(buffer)进行了加速处理













回复 支持 1 反对 0

使用道具 举报

11

主题

45

帖子

232

积分

中级会员

Rank: 3Rank: 3

积分
232
发表于 2017-4-2 08:18:38 | 显示全部楼层
x2yline 发表于 2017-3-27 19:45
077-x2yline
(1)处理fastq文件的python3代码如下:[mw_shl_code=python,true]# encoding=utf ...

你的代码写的很好、在buffer读取时剩余字符处理方式比我视频中的要简洁,代码复用性也不错,有一点需要注意一下,用buffer去读文件需要用大文件测试一下,读完再split速度并不一定有提高,这也是我没有用的原因,buffer更适合读到内存后不做其他处理的情况,希望你做个测试然后回复我,我也没有对这个做过测试的,文件用2G以上的
回复 支持 1 反对 0

使用道具 举报

29

主题

130

帖子

1136

积分

金牌会员

Rank: 6Rank: 6

积分
1136
发表于 2017-3-27 21:20:01 | 显示全部楼层

034-perl-shell

输入文件:longreads.fq
输出文件:longreads.fasta, longread_rm.fastq

代码如下:
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;

my $infile = "longreads.fq";
open my $fh, $infile or die "Cannot open $infile";
my ($total_len,$total_num,$total_GC,$max_num,$min_num,@array);
my $fasta = "longreads.fasta";
my $rm_fastq = "longread_rm.fastq";
open my $fh_out1, ">$fasta" or die "Cannot write to $fasta";
open my $fh_out2, ">$rm_fastq" or die "Cannot write to $rm_fastq";

my $trim_len = 5;
while(<$fh>){
	chomp;
	my $header = $_;
	chomp(my $seq = <$fh>);
	chomp(my $tmp = <$fh>);
	chomp(my $qual = <$fh>);
	my $len = length($seq);
	$total_len += $len;				#统计序列总长度以及碱基个数
	$total_num++;
	push @array, $len;
	my $GC = $seq =~ tr/GCgc/GCgc/;   
	$total_GC += $GC;				#统计GC含量
	print $fh_out1 ">$header\n$seq\n";  #FASTQ转变为FASTA
	my $rm_seq = substr($seq,$trim_len,$len-2*$trim_len);		#切去前后个5bp长度的碱基
	my $rm_qual = substr($qual,$trim_len,$len-2*$trim_len);
	print $fh_out2 "$header\n$rm_seq\n$tmp\n$rm_qual\n";
}
close $fh;
close $fh_out1;
close $fh_out2;

my $per_GC = sprintf "%2.2f", $total_GC/$total_len;
print "The number of base is $total_len\tthe content of GC is $per_GC\n";

my @seq_length = sort @array;
print "The total length of the reads is $total_len\n";
print "The min length of the reads is $seq_length[0]\n";
print "The min length of the reads is $seq_length[$#seq_length]\n";


回复 支持 1 反对 0

使用道具 举报

11

主题

45

帖子

232

积分

中级会员

Rank: 3Rank: 3

积分
232
发表于 2017-3-18 00:22:58 | 显示全部楼层
回复 支持 反对

使用道具 举报

0

主题

2

帖子

183

积分

注册会员

Rank: 2

积分
183
发表于 2017-3-27 14:14:43 | 显示全部楼层
本帖最后由 小清新范的吊丝 于 2017-3-27 14:18 编辑

Python - 294
解题思路:
对FASTQ的操作
    1、对文件进行循环读取,判断文件是否读取完毕。
    2、利用切片的方法,对5,3两端的数据进行截取
    3、每个循环都记下序列长度,最终统计得到总序列长度
    4、取fastq文件第一行和第二行得到fasta文件
    5、在循环中统计每个序列GC的含量,最后得到总的GC个数,从而得到GC%
[Python] 纯文本查看 复制代码
# -*- coding: utf-8 -*-
# @Time    : 2017/3/27 10:02
# @Author  : liang
# @File    : test1.0.py
# @Software: PyCharm
import sys
from optparse import OptionParser

args = sys.argv

def cut_fastq(filename,cut_5 = 0,cut_3 = 0):
    file_cut = open('read.fq', 'w')
    with open(filename) as f:
        while True:
            line_1 = f.readline().rstrip()
            if not line_1:
                break
            line_2 = f.readline().rstrip()[cut_5:][:-cut_3]
            line_3 = f.readline().rstrip()
            line_4 = f.readline().rstrip()[cut_5:][:-cut_3]

            file_cut.write(line_1 + '\n' + line_2 + '\n' + line_3 + '\n' + line_4 + '\n')

    file_cut.close()
def length_distribution(filename):
    dis = {}
    with open(filename) as f:
        while True:
            line_1 = f.readline().rstrip()
            if not line_1:
                break
            line_2 = f.readline().rstrip()
            line_3 = f.readline().rstrip()
            line_4 = f.readline().rstrip()
            length = len(line_2)
            if length not in dis:
                dis[length] = 0
            dis[length] += 1
    print dis

def fq_2_fa(filename):
    file_fasta = open('read.fa', 'w')
    with open(filename) as f:
        while True:
            line_1 = f.readline().rstrip()
            if not line_1:
                break
            line_2 = f.readline().rstrip()
            line_3 = f.readline().rstrip()
            line_4 = f.readline().rstrip()
            file_fasta.write('>'+line_1 + '\n' + line_2 + '\n')
    file_fasta.close()

def len_fq(filename):
    num = 0
    with open(filename) as f:
        while True:
            line_1 = f.readline().rstrip()
            if not line_1:
                break
            line_2 = f.readline().rstrip().upper()
            line_3 = f.readline().rstrip()
            line_4 = f.readline().rstrip()
            num = num + len(line_2)
    print num
    return num


def count_gc(filename):
    gc = 0
    with open(filename) as f:
        while True:
            line_1 = f.readline().rstrip()
            if not line_1:
                break
            line_2 = f.readline().rstrip().upper()
            line_3 = f.readline().rstrip()
            line_4 = f.readline().rstrip()
            gc = line_2.count('G') + line_2.count('C')
    SUM = len_fq(filename)
    print gc * 1.0 / SUM
    return gc * 1.0 / SUM


def main(args):
    parser = OptionParser()
    parser.add_option("-f", "--fastq", dest="filename",
                      help="fastq filename", metavar="FILE")

    parser.add_option("--cut", "--cut-fastq", dest="cut",
                      help="cut fastq file", action="store_true",
                      default=False)
    parser.add_option("-5", "--cut-5", dest="cut_5", type=int,
                      help="cut fastq N(N>0) bases from 5'", metavar="INT", default=0)
    parser.add_option("-3", "--cut-3", dest="cut_3", type=int,
                      help="cut fastq N(N>0) bases from 3'", metavar="INT", default=0)

    parser.add_option("--len-dis", dest="len_dis",
                      help="sequence length distribution", action="store_true",
                      default=False)

    parser.add_option("--fasta", "--fastq2fasta", dest="fasta",
                      help="convert fastq to fasta", action="store_true",
                      default=False)

    parser.add_option("--len", "--fastq-length", dest="fq_len",
                      help="total reads number", action="store_true",
                      default=False)

    parser.add_option("--gc", "--count-gc", dest="gc",
                      help="print GC%", action="store_true",
                      default=False)

    (options, args) = parser.parse_args()
    if not options.filename:
        parser.print_help()
    filename = options.filename

    if options.cut:
        cut_5 = options.cut_5
        cut_3 = options.cut_3
        cut_fastq(filename, cut_5, cut_3)

    if options.len_dis:
        length_distribution(filename)

    if options.fq_len:
        len_fq(filename)

    if options.fasta:
        fq_2_fa(filename)

    if options.gc:
        count_gc(filename)


if __name__ == '__main__':
    main(args)





小记:按照视频中的要求,我将剪切5,3两端序列之后的结果输出到文件 read.fq 文件中。将fq转换成fa文件存储到 read.fa 文件中。我对参数调用设定这个函数不熟悉,并没有掌握,需要加深学习。时间可能没有安排好,我想用这周时间对 optparse 模块加深学习。希望下次能完善好,能够使自己的代码重复利用起来。

对FASTA的操作
    1.由上面做的,我对fastq文件操作,产生一个fasta文件。用于这一步的编程。
    2.循环读取序列,并对序列取互补序列,并将产生的文件存放在 complement_reads.fa 文件中。考虑到有的序列是N,所以如果出现N,那么不操作。
    3.利用python中切片的方法直接将每个序列翻转顺序,得到反向序列。并将得到的反向序列存放在 reverse_reads.fa 文件中
    4.在每个循环中将T转换成U,得到RNA序列。并将RNA序列存放在 RNA.fa 文件中
    5.使用python中 readline(),确定size大小,从而输出制定长度的序列。
    6. 将每个ID和其序列对应为字典格式,这样可以输出制定ID的序列了。
    7.产生一个随机偶数,用来随机输出一个序列。
[Python] 纯文本查看 复制代码
# -*- coding: utf-8 -*-
# @Time    : 2017/3/27 11:35
# @Author  : liang
# @File    : test1.1.py
# @Software: PyCharm
import sys
import random

args = sys.argv

def complement_reads(filename):
    """
    取互补序列,并把生成的互补序列放在 complement_reads.fa 文件中
    """
    f1 = open('complement_reads.fa','w')
    with open(filename) as f:
        while True:
            line_1 = f.readline().rstrip()
            if not line_1:
                break
            line_2 = f.readline().rstrip().upper()
            for i in range(len(line_2)):
                if line_2[i] == 'A':
                    line_2[i] = 'T'
                elif line_2[i] == 'T':
                    line_2[i] = 'A'
                elif line_2[i] == 'C':
                    line_2[i] = 'G'
                elif line_2[i] == 'G':
                    line_2[i] = 'C'
                else:line_2[i] = 'N'
            f1.write(line_1+'\n'+line_2+'\n')
    f1.close()

def reverse_reads(filename):
    """
    取反向序列,并把生成的互补序列放在 reverse_reads.fa 文件中
    """
    f1 = open('reverse_reads.fa', 'w')
    with open(filename) as f:
        while True:
            line_1 = f.readline().rstrip()
            if not line_1:
                break
            line_2 = f.readline().rstrip()
            line_2 = line_2[::-1]
            f1.write(line_1 + '\n' + line_2 + '\n')
    f1.close()

def DNA_2_RNA(filename):
    """
    DNA序列转换成RNA序列,并放在RNA.fa文件中
    """
    f1 = open('RNA.fa', 'w')
    with open(filename) as f:
        while True:
            line_1 = f.readline().rstrip()
            if not line_1:
                break
            line_2 = f.readline().rstrip()
            for i in range(len(line_2)):
                if line_2[i] == 'T':
                    line_2[i] = 'U'
            f1.write(line_1 + '\n' + line_2 + '\n')
    f1.close()

def output_len(filename,length):
    """
    每行指定长度输出序列
    length 表示输出的长度
    """
    with open(filename) as f:
        while True:
            line_1 = f.readline().rstrip()
            if not line_1:
                break
            line_2 = f.readline(size=length).rstrip()
            print line_1
            print line_2

def ID_reads(filename,id_name):
    """
    提取指定文件指定ID的序列
    """
    dis = {}
    with open(filename) as f:
        while True:
            line_1 = f.readline().rstrip()
            if not line_1:
                break
            line_2 = f.readline().rstrip()
            dis[line_1[1:]] = line_2

    print dis[id_name]

def output_reads(filename):
    """
    随机抽取一个序列
    """
    with open(filename) as f:
        output_read = f.readlines()
        print output_read[random.randint(0,len(output_read)-1,2)]


小记: 由于我对optparse这个模块掌握并不明白,所以下面这个我只是把实现功能的函数写出来了,并没有整合后面的。这个模块我真的需要加强学习。 对FASTA文件操作中的 大小写字母形式输出我并不是很明白是什么意思,按照序列长度/名称排序 这个我并不会。所以这两个功能还没有做。
总之第一次做作业,时间没有安排好,所以做的很匆忙,通过看视频我也学到了很多别人编程的优点和长处,这对我非常非常有帮助。我觉得我要学习的还很多,这次没有安排好时间,下次一定安排好时间,把作业做好。
回复 支持 反对

使用道具 举报

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
发表于 2017-3-27 16:13:20 | 显示全部楼层

python刚开始学两天,比较费劲,一点一点的提交吧。当做每日笔记了。每天学一点。
加油
[AppleScript] 纯文本查看 复制代码
import sys#引入sys,为了获取参数
args=sys.argv#用于接下来获取命令行参数,我觉得这样赋值为了方便
filename=args[1]#第一个参数为文件名
cut_5=int(args[2])#由于参数为字符串所以改为int类型,否则会出错
cut_3=int(args[3])#参数获取完毕
  #读取文件
with open(filename) as FIN:#由于fastq文件是四行为一组的,不能用for  in 读取
    while TRUE:#用while一次读取四行
        line_1=FIN.readline().strip('\n')#读取4次
        if not(line_1 and FIN):#检查文件结束,防止出错
            break
        line_2=FIN.readline().strip('\n')
        line_3=FIN.readline().strip('\n')
        line_4=FIN.readline().strip('\n')
#读取序列,第2行是序列,第4行是质量,字符数相等
        print line_1
        print line_2 [cut_5:][:-cut_3]#3'段从5'段开始。由于5'段是反向,故用“-”号
        print line_3
        print line_4 [cut_5:][:-cut_3]

当每日笔记吧。今天还有别的事。明天再加。有错误请指出
回复 支持 反对

使用道具 举报

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
发表于 2017-3-27 18:19:13 | 显示全部楼层
[AppleScript] 纯文本查看 复制代码
import sys#引入sys,为了获取参数

args=sys.argv#用于接下来获取命令行参数,我觉得这样赋值为了方便

filename=args[1]#第一个参数为文件名

cut_5=int(args[2])#由于参数为字符串所以改为int类型,否则会出错

cut_3=int(args[3])#参数获取完毕

  #读取文件

with open(filename) as FIN:#由于fastq文件是四行为一组的,不能用for  in 读取

    while TRUE:#用while一次读取四行

        line_1=FIN.readline().strip('\n')#读取4次

        if not(line_1 and FIN):#检查文件结束,防止出错

            break

        line_2=FIN.readline().strip('\n')

        line_3=FIN.readline().strip('\n')

        line_4=FIN.readline().strip('\n')

#读取序列,第2行是序列,第4行是质量,字符数相等

        print line_1

        print line_2 [cut_5:][:-cut_3]#3'段从5'段开始。由于5'段是反向,故用“-”号

        print line_3

        print line_4 [cut_5:][:-cut_3]
回复 支持 反对

使用道具 举报

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
发表于 2017-3-27 18:32:18 | 显示全部楼层
学习最快乐 发表于 2017-3-27 16:13
python刚开始学两天,比较费劲,一点一点的提交吧。当做每日笔记了。每天学一点。
加油
[mw_shl_code=appl ...

希望大神原谅,一定好好学习。另外我的电脑卡了,然后就发了乱七八糟,
回复 支持 反对

使用道具 举报

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
发表于 2017-3-27 19:07:29 | 显示全部楼层
我懂的就只有这些了。以后会努力的。
[AppleScript] 纯文本查看 复制代码
#长度分布统计
import sys#引入sys,为了获取参数
args=sys.argv#用于接下来获取命令行参数
filename=args[1]#第一个参数为文件名
cut_5=int(args[2])#由于参数为字符串所以改为int类型,否则会出错
cut_3=int(args[3])#参数获取完毕

dis ={}#设置字典
  #读取文件
with open(filename) as FIN:#由于fastq文件是四行为一组的,不能用for  in 读取
    while TRUE:#用while一次读取四行
        line_1=FIN.readline().strip('\n')#读取4次
        if not(line_1 and FIN):#检查文件结束,防止出错
            break
        line_2=FIN.readline().strip('\n')
        line_3=FIN.readline().strip('\n')
        line_4=FIN.readline().strip('\n')

        length=len(line_2)
        if not length in dis:#判断字典内有没有
            dis[length]=0 #如果没有则赋值为0
            dis[length]+=1#有加1


for length ,num in dis.items():
                print("%s\t%s"%(length,num))
#fastq转换为fasta
import sys#引入sys,为了获取参数
args=sys.argv#用于接下来获取命令行参数
filename=args[1]#第一个参数为文件名
cut_5=int(args[2])#由于参数为字符串所以改为int类型,否则会出错
cut_3=int(args[3])#参数获取完毕
  #读取文件
dis={}
with open(filename) as FIN:#由于fastq文件是四行为一组的,不能用for  in 读取
    while TRUE:#用while一次读取四行
        line_1=FIN.readline().strip('\n')#读取4次
        if not(line_1 and FIN):#检查文件结束,防止出错
            break
        line_2=FIN.readline().strip('\n')
        line_3=FIN.readline().strip('\n')
        line_4=FIN.readline().strip('\n')
#只在第一行加>,第二行读出就好
        print'>'+line_1
        print line_2
        
#计算长度以及GC含量
import sys#引入sys,为了获取参数
args=sys.argv#用于接下来获取命令行参数
filename=args[1]#第一个参数为文件名
cut_5=int(args[2])#由于参数为字符串所以改为int类型,否则会出错
cut_3=int(args[3])#参数获取完毕
  #读取文件
NUM=0#设置初始值,下同
GC=0
with open(filename) as FIN:#由于fastq文件是四行为一组的,不能用for  in 读取
    while TRUE:#用while一次读取四行
        line_1=FIN.readline().strip('\n')#读取4次
        if not(line_1 and FIN):#检查文件结束,防止出错
            break
        line_2=FIN.readline().strip('\n').upper()
        line_3=FIN.readline().strip('\n')
        line_4=FIN.readline().strip('\n')

        NUM+=len(line_2)#因为求的是序列长度,所以只用计算第二行即可
        GC+=line_2.count('G')+line_2.count('C')#GC含量,同样也是只计算第二行

        
print NUM
print GC*1.0/NUM  #*1.0是为了有小数,否则就是个整型数    

最后如果是版本3,print后加()
知道菜的不能再菜了,但是希望自己不会掉队的心,还是很像交作业的,虽然犯很多蠢蠢的问题。希望见谅。我会好好学python的。
回复 支持 反对

使用道具 举报

0

主题

1

帖子

29

积分

新手上路

Rank: 1

积分
29
发表于 2017-3-27 21:14:29 | 显示全部楼层
本帖最后由 张起灵 于 2017-3-27 21:27 编辑

[Python] 纯文本查看 复制代码
#encoding=utf-8

"""
简介: fastq文件数据处理工具(因为觉得Biopython很好用,便用biopython完成了作业)
作者: Xiaoxi Zhang
时间:2017年3月25日
"""
import sys
import  Bio
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from optparse import OptionParser
import csv

args=sys.argv

def cut_fastq(fq_file, fq_cut_file,cut_5, cut_3):
    """
    从fastq文件的5'端 3'端截掉几个碱基并存入本地fastq文件,因为更多的时候是截取的序列拿来做后续分析,
    所以就直接存起来而不输出到屏幕上了。另外本来想用SeqIO读取fastq文件来完成这一题的,但是奈何在
    biopython里面是用嵌套的字典来处理质量值的,所以处理起来好麻烦,还不如老老实实一行一行读取来得快
    :param fq_file:需要截取的fastq文件名
    :param fq_cut_file:存入本地的截取后的fastq文件名
    :param cut_5:5'端截取的碱基数,大于或等于0
    :param cut_3:3'端截取的碱基数,大于或等于0
    :return:
    """
    if cut_3<0:
        raise ValueError('cut_3: %s < 0'% cut_3)
    if cut_5<0:
        raise ValueError('cut_5: %s < 0'% cut_5)
    with open (fq_file) as FIN:
        while FIN:
            line_1=FIN.readline()
            if not (line_1 and FIN):
                break
            line_2=FIN.readline()
            line_3=FIN.readline()
            line_4=FIN.readline()
            fq_file=open(fq_cut_file,"a") #使用追加模式,不然后面的会覆盖前面写入的内容
            fq_file.write(line_1)
            if cut_5>=0 and cut_3>0:  #当有时候只需要截取一端的时候
                fq_file.write(line_2[cut_5:-cut_3])
            else:
                fq_file.write(line_2[cut_5:])
            fq_file.write(line_3)
            fq_file.write(line_4[cut_5:-cut_3])

def length_distribution(fq_file,csv_file):
    """
    统计fastq文件的长度分布输出到屏幕并存入到本地csv文件方便后续画长度分布图
    :param fq_file: fastq文件名
    :param out_file: 存入本地的csv文件名
    :return:
    """
    length_list=[len(seq) for seq in SeqIO.parse(fq_file,"fastq")]
    length_set=set(length_list) #length_set是另外一个列表,里面的内容是length_list里面的无重复项
    csvfile=csv.writer(open(csv_file,"w",newline=''))  ##创建并以写入模式打开一个csv文件
    for length in length_set:
        print(length,length_list.count(length))
        csvfile.writerow([length,length_list.count(length)])
    return length_list

def fq_2_fa(fq_file,fa_file):
    """
    将fastq文件转换为fasta文件
    :param fq_file: fastq文件名
    :param fa_file: fasta文件名
    :return:
    """
    #SeqIO.write(list(SeqIO.parse(fq_file,"fastq")),fa_file,"fasta")  这样写可以 但是远没有现成的转化函数好用
    #转化函数可以转化很多的序列格式,像gbk,fastq,fasta等等一行代码实现互相转化
    SeqIO.convert(fq_file,"fastq",fa_file,"fasta")

def fq_length(fq_file):
    """
    fastq文件总碱基数统计
    :param fq_file: fastq文件名
    :return:
    """
    length=sum([len(seq) for seq in SeqIO.parse(fq_file,"fastq")])
    #其实这一步可以调用length_distribution()函数的返回值
    print(length)
    return length

def count_gc(fq_file):
    """
    统计fastq文件的GC总碱基数
    :param fq_file: fastq文件名
    :return:
    """
    gc_sum=0
    for fq_seq in SeqIO.parse(fq_file,"fastq"):
        gc_sum+=fq_seq.seq.upper().count("G")+fq_seq.seq.upper().count("C")
    print("GC percentage: %s"% (gc_sum*1.0/fq_length(fq_file)))

#作为一名菜鸡,之前完全不知道有这个OptionParser模块,之前还一直好奇
# 服务器上面那些帮助文档是怎么写出来,原来是这样来的,好玩!高级!
def main(args):
    usage = "usage: %prog [options] arg1 arg2..."
    parser= OptionParser(usage)
    parser.add_option("-q", "--fastq", dest="fq_file", help="fastq filename",metavar="FILE")
    parser.add_option("-C", "--cut-fastq", dest="cut", help="cut fastq file",
                      action="store_true", default=False, metavar="FILE")
    parser.add_option("-c", "--output-cut-fastq-file", dest="fq_cut_file", help="output fastq filename when cut bases from 5' or 3'"
                      ,metavar="FILE")
    parser.add_option("-5", "--cut-5", dest="cut_5", type=int, help="cut fastq N(N>0) bases from 5 end",
                      metavar="INT", default=0)
    parser.add_option("-3", "--cut-3", dest="cut_3", type=int, help="cut fastq N(N>0) bases from 3 end",
                      metavar="INT", default=0)

    parser.add_option("-D", "--length-dis", dest="len_dis", help="analysis fastq length distribution",
                      action="store_true",default=False)
    parser.add_option("-v", "--csv", dest="csv_file", help="output csv filename when analysis reads length distribution"
                      ,metavar="FILE")

    parser.add_option("-2", "--fastq2fasta", dest="fq_fa", help="convert fastq to fasta",
                      action="store_true",default=False)
    parser.add_option("-a", "--fasta", dest="fa_file", help="output fasta filename when convert fastq to fasta",metavar="FILE")

    parser.add_option("-L", "--fastq-length", dest="fq_len", help="total reads number",
                      action="store_true",default=False)

    parser.add_option("-G", "--count-gc", dest="gc",help="print GC%",action="store_true",default=False)

    (options,args)=parser.parse_args()  #解析args命令行参数
    if not options.fq_file:
        parser.print_help()
    fq_file=options.fq_file
    fa_file=options.fa_file
    csv_file=options.csv_file
    fq_cut_file=options.fq_cut_file

    if options.cut:
        cut_5=options.cut_5
        cut_3=options.cut_3
        cut_fastq(fq_file,fq_cut_file,cut_5,cut_3)

    if options.len_dis:
        length_distribution(fq_file,csv_file)

    if options.fq_fa:
        fq_2_fa(fq_file,fa_file)

    if options.fq_len:
        fq_length(fq_file)

    if options.gc:
        count_gc(fq_file)
if __name__=='__main__':
    main(args)


回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2018-5-21 05:45 , Processed in 0.108801 second(s), 28 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.