搜索
楼主: Jimmy

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

[复制链接]

29

主题

131

帖子

1208

积分

金牌会员

Rank: 6Rank: 6

积分
1208
发表于 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

使用道具 举报

0

主题

4

帖子

89

积分

注册会员

Rank: 2

积分
89
QQ
发表于 2017-3-28 20:49:13 | 显示全部楼层
本帖最后由 xzmcxjb 于 2017-3-31 11:49 编辑

徐**连二ICU作业,
1.第一周第四题统计碱基个数及GC含量python编码及结果

import sys#调用系统信息相关模块

args=sys.argv#命名sys模块参数

filename=args[1]#定义文件名
cut_5=int(args[2])#将args[2]文件转换为数值型整型,命名为cut_5
cut_3=int(args[3])#将args[3]文件转换为数值型整型,命名为cut_3


NUM=0#将NUM值清零
GC=0#将GC值清零
with open(filename) as FIN:#调用FIN文件
    while True: #无限循环模式
        line_1 = FIN.readline().strip('\n')#读取第一行去除换行
        if not(line_1 and FIN):
            break                                     #如果不是FIN文件第一行则中断
        line_2 = FIN.readline().strip('\n').upper()#读取第二行去除换行,采用大写字符
        line_3 = FIN.readline().strip('\n')#读取第三行去除换行
        line_4 = FIN.readline().strip('\n')#读取第四行去除换行

        NUM+=len(line_2)                   #第二行碱基综述命名为NUM
        GC+=line_2.count('G')+line_2.count('C')#分别计算第二行G与C的数量,命名为GC
print (NUM)                                      #输出NUM
print (GC*1.0/NUM)                          #输出GC数与NUM的比值
结果:
C:\Users\Administrator\AppData\Local\Programs\Python\Python36\python.exe C:/Users/Administrator/PycharmProjects/untitled1/biotrainee0.py reads_1.fq 3 5
1088399
0.4869381541144378

Process finished with exit code 0
第二单元第四题大小写字母形式输出
import sys#调用sys模块
args=sys.argv#命名sys模块参数
filename=args[1]#命名filename
out_filename=1#先在python程序文件夹建立一个文本文件1.txt,命名为out_filename
with open(filename) as F_in:#指定filename为输入文件
    with open(out_filename, 'w') as F_out:#指定out_filename为输出文件
        for line in F_in:
            line = line.strip('\n')#将F_in里的每行去掉换行
            if not line.startswith('>'):#如果序列不是FASTA格式
                if capital:
                    seq = line.upper()#如果行内有大写字符,一律大写
                else:
                    seq=line.lower()#否则小写
            else:
                seq = line#如果序列是FASTA格式
            F_out.write('%s\n' % (seq))#直接输出序列

第二周第五题实在不懂,容慢慢领会。



回复 支持 反对

使用道具 举报

0

主题

2

帖子

123

积分

注册会员

Rank: 2

积分
123
发表于 2017-3-30 15:48:11 | 显示全部楼层
想法:把原文件的fastq写入list进行操作;python画图不是很懂,借用了!!!077-x2yline!!!的可视化代码!(肉眼不可见的改动);用了另一种argparse作为函数,感觉很好理解;
代码如下:
[Python] 纯文本查看 复制代码
#!/usr/bin/env python
# -*- coding: utf-8 -*-

# ---------------------------------------------------------------
# Several tools for sequences 
# ---------------------------------------------------------------

"""
seqtools.py
============
Please type "./seqtools.py -h" for usage help
    
Author:
    Chunyu

Description:
    Several tools for sequences
Reurirements:
    Python packages: argparse,matplotlib,scipy,os,numpy
"""

# ---------------------------------------------------------------
# import
# ---------------------------------------------------------------

import argparse
from matplotlib import pyplot as plt
from scipy import stats
import os
import numpy as np
# ---------------------------------------------------------------
# function definition
# ---------------------------------------------------------------

def parse_args():
    """ master argument parser """
    
    parser = argparse.ArgumentParser(
        description = "",
        #epilog="",
        #formatter_class=argparse.RawDescriptionHelpFormatter,
        )
    parser.add_argument(
        '-i', '--input_file_fq',
        type=str,
        required=True,
        help="""
        Transmit file inputted to the script.
        """,
        )
    parser.add_argument(
        '-o', '--output_file',
        type=str,
        #required=True,
        default="default.output",
        help="""
        Define name of file outputted.
        """,
        )
    parser.add_argument(
        '-cut5', '--cut_5',
        type=int,
        #required=True,
        default=0,
        help="""
        Define the number of 5'cut.
        """,
        )
    parser.add_argument(
        '-cut3', '--cut_3',
        type=int,
        #required=True,
        default=0,
        help="""
        Define the number of 3'cut.
        """,
        )
    args = parser.parse_args()
    return args
def getfile(input_file_fastq):
    count=0
    file_list_fa=[]
    file_list_fq=[]
    with open (input_file_fastq,'r') as input_file:
        while input_file:
            line_1=input_file.readline()
            if not line_1:
                break
            line_1.strip()
            line_2=input_file.readline()
            line_2=file_list_fa.append(line_2.strip().upper())
            line_3=input_file.readline().strip()
            line_4=input_file.readline()
            line_4=file_list_fq.append(line_4.strip().upper())
            #print file_list_fa
    return file_list_fa,file_list_fq

def cut_53(file_list_fa,cut_5,cut_3):
    file_list_cut53=[]
    for i in file_list_fa:
        file_list_cut53.append(i[cut_5:][:-cut_3])
    return file_list_cut53   

def length_distribution(file_list_fa): 
    length=[]
    for i in file_list_fa:
        length.append(len(i))
        
    return length

def gc_content(file_list_fa):
    GC_all_number=0
    N_number=0
    length=0
    GC_distribution=[]
    for i in file_list_fa:
        length=length+len(i)
    for i in file_list_fa:
        GC_tmp=i.count('C')+i.count('G')
        GC_number=GC_all_number+GC_tmp
        N_number=N_number+i.count('N')
        GC_distribution.append(round((GC_tmp*1.0)/(len(i)*1.0),3))
        #print GC_tmp
        #print round((GC_tmp*1.0)/(len(i)*1.0),3)
    return GC_number/length,N_number/length,GC_distribution

def complement(file_list_fa):
    ACTG_dic={'A':'T','T':'A','C':'G','G':'C','N':'N'}
    TGAC=[]
    for i in file_list_fa:
        for j in i:
            TGAC.append(ACTG_dic[j])
    return TGAC

def reverse(file_list_fa):
    ACTG=[]
    for i in file_list_fa:
        ACTG.append(i)
    GCAT=''.join(ACTG)[::-1]
    return GCAT

def DNA2RNA(file_list_fa):
    file_list_reverse=[]
    for i in file_list_fa:
        try:
            file_list_reverse.append(i.replace('T','U'))
        except:
            pass
    return file_list_reverse

def length_rank(file_list_fa):
    file_list_fa.sort(key=lambda x:len(x))
    return file_list_fa
        
            
def density_plot(data_list,xlabel,title):
    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([min(data_list)-data_range/10,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(min(data_list)-data_range/10,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))))                                       
                

    
# ---------------------------------------------------------------
# main function
# ---------------------------------------------------------------
def main():
    """ main """

    args = parse_args()
    input_file_fastq=args.input_file_fq
    output_file=args.output_file
    with open(output_file,'w+') as file_output:
        #put original file into a list as fasta
        fasta_list,fastq_list=getfile(input_file_fastq)
        cut_5=args.cut_5
        cut_3=args.cut_3
        file_output.writelines("This is the result after cutting 5'& 3'"+'\n')
        file_output.writelines(cut_53(fasta_list,cut_5,cut_3))
        file_output.writelines('\n')
        file_output.writelines("This is the result of complement"+'\n')
        file_output.writelines(complement(fasta_list))
        file_output.writelines('\n')
        length=length_distribution(fasta_list)
        density_plot(length,'length','length_distribution')
        GC_percent,N_percent,GC_distribution=gc_content(fasta_list)
        file_output.writelines("This is the result of GC_distribution"+'\n')
        file_output.writelines(str(GC_distribution)[1:-1])
        density_plot(GC_distribution,'gc_number','gc_distribution')
        file_output.writelines('\n')
        file_output.writelines("This is the result of reverse"+'\n')
        file_output.writelines(reverse(fasta_list))
        file_output.writelines('\n')
        file_output.writelines("This is the result of DNA2RNA"+'\n')
        file_output.writelines(DNA2RNA(fasta_list))
        file_output.writelines('\n')
        file_output.writelines("This is the result of length_rank"+'\n')
        for i in length_rank(fasta_list):
            file_output.writelines(i)
            file_output.writelines('\n')

if __name__ == "__main__":
    main()

回复 支持 反对

使用道具 举报

0

主题

5

帖子

105

积分

注册会员

Rank: 2

积分
105
发表于 2017-3-31 17:21:38 | 显示全部楼层
本帖最后由 王土土哈 于 2017-3-31 17:27 编辑

2nd fasta文件操作
[Python] 纯文本查看 复制代码
# -*- coding: utf-8 -*-
"""
Created on Fri Mar 31 15:06:18 2017

@author: moon
"""
import os
os.chdir("E:\\work\\biopython\\biotrainee\\lesson0")
from Bio.SeqRecord import SeqRecord
from Bio import SeqIO

#大小写字母形式输出
def lower_upper_switch(filename, to_upper = 0, to_lower = 0):
    if to_upper == 1:
        if to_lower == 1:
            print("arguments wrong! We can't not print 2 kinds of fasta.")
        else:
            records = (rec.upper() for rec in SeqIO.parse(filename, "fasta"))
            SeqIO.write(records, "upper.fas", "fasta")
    elif to_lower == 1:
        records = (rec.lower() for rec in SeqIO.parse(filename, "fasta"))
        SeqIO.write(records, "lower.fas", "fasta")
    else:
        print("No switch of fasta type.")
        
lower_upper_switch("ath_miRNA.fa")

#DNA to RNA
def make_RNA_record(nuc_record):
    """Returns a new SeqRecord with the transcribed sequence."""
    return SeqRecord(seq = nuc_record.seq.transcribe(), \
                     id = "trans_" + nuc_record.id, \
                     description = "transcribe from DNA")
def DNA2RNA(filename, transcribe = "N"):
    if transcribe == "Y":
        records = (make_RNA_record(nuc_rec) for nuc_rec in SeqIO.parse(filename, "fasta"))
        SeqIO.write(records, "transcribe.fasta", "fasta")
    else:
        print("No change.")
    
DNA2RNA("ath_miRNA.fa",transcribe = "Y")

#取互补序列
def make_complement_record(nuc_record):
    """Returns a new SeqRecord with the complement sequence."""
    return SeqRecord(seq = nuc_record.seq.complement(), \
                     id = "complement_" + nuc_record.id, \
                     description = "complement from DNA")
def DNA_complement(filename, complement = "N"):
    if complement == "Y":
        records = (make_complement_record(nuc_rec) for nuc_rec in SeqIO.parse(filename, "fasta"))
        SeqIO.write(records, "complement.fasta", "fasta")
    else:
        print("No change.")
    
DNA_complement("ath_miRNA.fa",complement = "Y")

#取反向序列
def DNA_reverse(filename, reverse = "N"):
    if reverse == "Y":
        records = (nuc_rec[::-1] for nuc_rec in SeqIO.parse(filename, "fasta"))
        SeqIO.write(records, "reverse.fasta", "fasta")
    else:
        print("No change.")
    
DNA_reverse("ath_miRNA.fa",reverse = "Y")





回复 支持 反对

使用道具 举报

11

主题

52

帖子

279

积分

中级会员

Rank: 3Rank: 3

积分
279
发表于 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

使用道具 举报

11

主题

52

帖子

279

积分

中级会员

Rank: 3Rank: 3

积分
279
发表于 2017-4-2 08:27:20 | 显示全部楼层
张起灵 发表于 2017-3-27 21:14
[mw_shl_code=python,true]#encoding=utf-8

"""

Bio.Seq的包确实不错、我建议你去看一下每个功能的源码来学习编程技巧,最后自己写一个类似的包给自己用或者对Bio包做扩展,例如把GC含量的功能加进去
回复 支持 反对

使用道具 举报

11

主题

52

帖子

279

积分

中级会员

Rank: 3Rank: 3

积分
279
发表于 2017-4-2 08:33:50 | 显示全部楼层
生信小菜鸟 发表于 2017-3-30 15:48
想法:把原文件的fastq写入list进行操作;python画图不是很懂,借用了!!!077-x2yline!!!的可视化代码 ...

知道lambda是怎么工作的么、可以查一下,python里很重要的一个特性
回复 支持 反对

使用道具 举报

2

主题

34

帖子

773

积分

高级会员

Rank: 4

积分
773
发表于 2017-4-2 22:50:29 | 显示全部楼层
本帖最后由 x2yline 于 2017-4-3 09:58 编辑
dongye 发表于 2017-4-2 08:18
你的代码写的很好、在buffer读取时剩余字符处理方式比我视频中的要简洁,代码复用性也不错,有一点需要注 ...

测试了一下,发现一个bug,在读取buffer后split('\n@')时有可能并不是每个read的分界(有可能是准确率那一行的首个字母是@),所以后来又改动了一下,用split('\n@SRR')进行分割。
测试的时候,发现大数据进行绘图耗时不是一般的长,所以我把可视化的代码给注释掉后进行了测试比较。
所用数据为ftp://ftp-trace.ncbi.nih.gov/sra ... 6632/SRR5396632.sra转化而成的SRR5396632.fastq,fastq文件大小为10.8G,运行速度:


项目 readline() buffer
GC含量1217.0491607189178 s 293.6417953968048 s
长度统计307.22693967819214 s243.5819320678711 s
fq2fa1147.9280161857605 s 466.1572198867798 s
cut 4 5   1558.3281314373016 s 672.0114369392395 s
ps:由于我的fq2fa程序把序列名中的@去掉了,所以生成的fasta文件比dongye老师的程序结果小一点点。

回复 支持 反对

使用道具 举报

29

主题

131

帖子

1208

积分

金牌会员

Rank: 6Rank: 6

积分
1208
发表于 2017-4-4 02:54:08 | 显示全部楼层
本帖最后由 anlan 于 2017-4-4 02:56 编辑

034 -perl+shell

FASTQ文件探究part2-Per:视频的那个抽取随机序列的方法之前试过,如果序列量上百万的话,个人电脑会很慢,之后用了下述方法,速度提升不少。
代码如下:
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;
my $in = "longreads.fasta";

###取互补序列###
open my $fh, $in or die "Cannot open $in";
open my $fh1, ">com_longreads.fasta" or die "Cannot write it";
while(<$fh>){
        chomp;
        if($_ =~ />/){
                print $fh1 "$_\n";
        }else{
                $_ =~ tr/ATCGatcg/TAGCtagc/;
                print $fh1 "$_\n";
        }
}
close $fh;
close $fh1;

###取反向序列###
open $fh, $in or die "Cannot open $in";
open my $fh2, ">reverse_longreads.fasta" or die "Cannot write it";
while(<$fh>){
        chomp;
        if($_ =~ />/){
                print $fh2 "$_\n";
        }else{
                my $tmp = reverse $_;
                print $fh2 "$tmp\n";
        }        
}
close $fh;
close $fh2;

###DNA to RNA###
open $fh, $in or die "Cannot open $in";
open my $fh3, ">rna_longreads.fasta" or die "Cannot write it";
while (<$fh>){
        chomp;
        chomp;
        if($_ =~ />/){
                print $fh3 "$_\n";
        }else{
                $_ =~ s/T/U/g;
                print $fh3 "$_\n";
        }        
}
close $fh;
close $fh3;

###每行指定长度输出###
open $fh, $in or die "Cannot open $in";
open my $fh4, ">longreads2.fasta" or die "Cannot write it";
while(<$fh>){
        chomp;
        my $len = length($_);
        if ($_ =~ />/){
                print $fh4 "$_\n";
        }else{
                for (my $i=0;$i<$len;$i+=70){
                        my $seq = substr($_,$i,70);
                        print $fh4 "$seq\n";
                }
        }
}
close $fh;
close $fh4;

###随机抽取序列###
open $fh, $in or die "Cannot open $in";
open my $fh5, ">sample_longreads.fasta" or die "Cannot write it";
my %hash;
my $header;
while(<$fh>){
        chomp;
        if ($_ =~ />/){
                $header = $_;
        }else{
                $hash{$header} = $_;
        }
}

my $count = keys %hash;
my @id = keys %hash;
my $k = int($count/10);
my @array = &random_num($count, $k);
map{
        print $fh5 "$id[$_-1]\n$hash{$id[$_-1]}\n";
}@array;


sub random_num {
        my ($k,$count) = @_;
        my $i;
        my @array = 1..$k;
        for ($i=0;$i<$count;$i++){
                my $index = int(rand($k-1));
                my $a = $array[$i];
                my $b = $array[$index];
                $array[$index] = $a;
                $array[$i] = $b;
        }
        my $j = 0;
        my @tmp;
        while($j<$count){
                push @tmp, $array[$j];
                 $j++;
        }
        return @tmp;
}
close $fh;
close $fh5;







回复 支持 反对

使用道具 举报

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

使用道具 举报

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

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-21 01:15 , Processed in 0.036194 second(s), 22 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.