搜索
查看: 4200|回复: 3

生信编程直播第六题:批量根据基因list来提取信息(shell)

[复制链接]

633

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2017-1-27 11:12:32 | 显示全部楼层 |阅读模式
如果只有一个基因list,那么从大的文件里面提取信息,很容易。比如:详细讲解如何抽取fasta文件里面指定序列名的序列 http://www.biotrainee.com/thread-150-1-1.html 

如果有多个多个gene list,就需要写循环啦,这个时候就用shell啦。
[Bash shell] 纯文本查看 复制代码
ls *txt |while read id
do 
cat $id ~/annotation/CHIPseq/mm10/ucsc.refseq.bed |perl -alne '{$h{$F[0]}=1;print if exists $h{$F[3]} }' >${id%%.*}.bed
done


请自行去UCSC下载refseq的bed,然后自己找几个gene list去这个refseq.bed里面提取自己的gene list的bed文件。



上一篇:全基因组分窗口统计GC含量
下一篇:【直播】我的基因组(七):从整体理解全基因组测序数据...
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

1

主题

15

帖子

180

积分

注册会员

Rank: 2

积分
180
发表于 2017-2-26 12:11:31 | 显示全部楼层
你的方法好巧妙 学习了
[Shell] 纯文本查看 复制代码
for i in *.txt                                          
cat $i |while read in
do
grep "$in" ref.bed >>$i.out
done
回复 支持 反对

使用道具 举报

0

主题

8

帖子

119

积分

注册会员

Rank: 2

积分
119
发表于 2017-6-30 10:07:13 | 显示全部楼层
这题也可以用Python来实现,代码复杂了点,但是可以用:
[Python] 纯文本查看 复制代码
import glob
import sys
args=sys.argv

#how to write the print result to the coordinated file



def readID(file):#store the ID of all files to a list
	idList=[]#this list should be defined in this sub function
	with open(file) as f:
		for line in f:
			line=line.strip('\n')
			idList.append(line)
	return idList		
def main(args):
	txtList=glob.glob(r'id*.txt')
	print(txtList)#to show the txt file the function find
	outfilename='out'#define the former part of the output file name
	for txtfile in txtList:

		#idList=[]
		acceptList=readID(txtfile)#accept the return value of readID
		print(acceptList)#show the list accepted from the sub function
		with open(args[1]) as F_in:
			with open("%s_%s"%(outfilename,txtfile),'w') as F_out:#two or more format name should be used as tuple
				for line in F_in:
					key=line.split('\t')[3]
					if key in acceptList:
						F_out.write(line.strip('\n')+'\n')#the line contains the '\n',so in fact we needn't addd '\n'
						print(line)
if __name__=="__main__":#the entry of the main
	main(args)
回复 支持 反对

使用道具 举报

2

主题

23

帖子

169

积分

注册会员

Rank: 2

积分
169
发表于 2019-6-10 00:10:55 | 显示全部楼层
boyxujingen 发表于 2017-6-30 10:07
这题也可以用Python来实现,代码复杂了点,但是可以用:
[mw_shl_code=python,true]import glob
import sys ...

你好:
能不能分享一下‘refseq.bed’文件,在UCSC没找到!谢谢
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-11-19 22:19 , Processed in 0.037254 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.