查看: 252|回复: 6

[Linux系统] 如何从fasta文件中提取大于一定长度的序列?

[复制链接]
  • TA的每日心情

    2017.8.15 10:27
  • 签到天数: 29 天

    连续签到: 1 天

    [LV.4]偶尔看看III

    帝王蝶

    Rank: 4

    主题
    9
    奥币
    514
    积分
    327
    注册时间
    2016.4.8
    在线时间
    31 小时

    发表于 2019.8.19 09:06:47 | 显示全部楼层 |阅读模式
    20奥币
    如何从fasta文件中提取大于一定长度的序列?或者哪个软件能做到这一点?比如10万个序列,长度不一,从中提取>300bp的所有序列。多谢!

    回复

    使用道具 举报

  • TA的每日心情
    吃饭
    2019.7.22 10:38
  • 签到天数: 450 天

    连续签到: 1 天

    [LV.9]以坛为家II

    版主

    Rank: 10Rank: 10Rank: 10

    主题
    21
    奥币
    4712
    积分
    1860
    注册时间
    2015.12.29
    在线时间
    431 小时

    突出贡献优秀版主热心会员


    发表于 2019.8.19 10:06:25 | 显示全部楼层
    回复

    使用道具 举报

  • TA的每日心情
    忙~
    昨天 09:27
  • 签到天数: 256 天

    连续签到: 7 天

    [LV.8]以坛为家I

    帝王蝶

    Rank: 4

    主题
    11
    奥币
    1062
    积分
    430
    注册时间
    2018.1.3
    在线时间
    171 小时

    热心会员活跃会员


    发表于 2019.8.19 17:19:51 | 显示全部楼层
    所有序列在一个大文件里的话,用这个脚本可以输出你想要的序列,但是我把序列变成一行了,内部换行的话重新处理也很简单的
    [Python] 纯文本查看 复制代码
    import sys
    import re
    
    dct,chR,seq={},'',[]
    f=open(sys.argv[1],'r')
    for i in f:
            if i.startswith(">"):
                    chR=i[1:-1]
                    seq=[]
                    dct[chR]=seq
            else:
                    seq.append(i)
    
    for i in dct:
            seq2=re.sub("\n","","".join(dct[i]))
            if len(seq2)>300:
                    print(">"+i+"\n"+seq2)
    
    回复

    使用道具 举报

  • TA的每日心情

    2017.8.15 10:27
  • 签到天数: 29 天

    连续签到: 1 天

    [LV.4]偶尔看看III

    帝王蝶

    Rank: 4

    主题
    9
    奥币
    514
    积分
    327
    注册时间
    2016.4.8
    在线时间
    31 小时

     楼主| 发表于 2019.8.22 08:49:50 | 显示全部楼层
    Ivan 发表于 2019.8.19 10:06
    https://bioinf.shenwei.me/seqkit/

    你好,这个软件好像没有我要的功能,他有个功能是split,但是作用是split sequences into multi parts with N sequences,这个size的意思可能还是指序列的数量,不是大小,不过感谢回复,这个软件有个translate的功能,很好用!
    回复

    使用道具 举报

  • TA的每日心情

    2017.8.15 10:27
  • 签到天数: 29 天

    连续签到: 1 天

    [LV.4]偶尔看看III

    帝王蝶

    Rank: 4

    主题
    9
    奥币
    514
    积分
    327
    注册时间
    2016.4.8
    在线时间
    31 小时

     楼主| 发表于 2019.8.22 08:51:16 | 显示全部楼层
    大裤衩 发表于 2019.8.19 17:19
    所有序列在一个大文件里的话,用这个脚本可以输出你想要的序列,但是我把序列变成一行了,内部换行的话重新 ...

    多谢回复,不知道有没有软件可以做到这一点?
    回复

    使用道具 举报

  • TA的每日心情
    吃饭
    2019.7.22 10:38
  • 签到天数: 450 天

    连续签到: 1 天

    [LV.9]以坛为家II

    版主

    Rank: 10Rank: 10Rank: 10

    主题
    21
    奥币
    4712
    积分
    1860
    注册时间
    2015.12.29
    在线时间
    431 小时

    突出贡献优秀版主热心会员


    发表于 2019.8.22 10:41:47 | 显示全部楼层
    张俊亚 发表于 2019.8.22 08:49
    你好,这个软件好像没有我要的功能,他有个功能是split,但是作用是split sequences into multi parts wi ...

    请认真看一下,是可以实现你的要求。
    先统计序列长度:seqkit faidx  yourfa
    然后提取序列长度大于300bp的:  cat yourfa.fai |awk '$2>300{print $1}' > large300
    最后提取序列: seqkit grep  -f large300   yourfa >  large300.fa

    会写程序的话,也很简单。
    回复

    使用道具 举报

  • TA的每日心情
    吃饭
    2019.7.22 10:38
  • 签到天数: 450 天

    连续签到: 1 天

    [LV.9]以坛为家II

    版主

    Rank: 10Rank: 10Rank: 10

    主题
    21
    奥币
    4712
    积分
    1860
    注册时间
    2015.12.29
    在线时间
    431 小时

    突出贡献优秀版主热心会员


    发表于 2019.8.22 10:44:31 | 显示全部楼层
    张俊亚 发表于 2019.8.22 08:49
    你好,这个软件好像没有我要的功能,他有个功能是split,但是作用是split sequences into multi parts wi ...

    你测试一下就知道了,seqkit很好用,多测试一下里面的参数,你会有意外的收货。
    回复

    使用道具 举报

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

    本版积分规则

    快速回复 返回顶部 返回列表