提取指定区间的转录本序列
2019-03-30 本文已影响3人
chaimol
用户需求：给定一个区间，找到该区间内所有的基因、转录本。
获得所有的转录本序列。
解决方案：
1.获取目标区段内所有的基因、转录本
cat merged.gtf|awk '$1 ~/8/ && $4>90193841 && $5<91589697 {print $0}' >tran_8
tran_8即是给出所有的在区间内的基因、转录本。
参数讲解：
awk $1匹配第一列 以8开头的，此处是第8染色体。 && &运算符，同时满足三个条件的，返回后面指定的列。
2.从基因组提取指定的序列。
cat tran_8 |awk '$2 ~/trans/ {print $8}|wc -l'
统计总共有多少个转录本。
实现方案：
- bedtools（常用工具）
- pysam (python包）
先用bedtools，据说速度飞快。
参考链接
说明一点：bedtools其中一个输入文件格式要求是tab分隔符的，我的解决方案是先用awk选好文件内容，之后导入本地的excel，另存为制表符格式。否则，bedtools会一直报错，说你的输入不够3列
我太喜欢awk,sed,grep这种神奇的命令，快捷干净，无废话。
1.先从所有的文件中找出你需要的区间内的所有转录本位置
cat merged.gtf|awk '$1 ~/8/ && $4>90193841 && $5<91589697 {print $0}' >tran_8
-
导入本地，转换为制表符格式。注意，尽量格式顺序如下
image.png
第4列是你要的命名的转录本名称。把第9列和第5列合并，即为第4列。第一列是chr8还是8根据你的基因组序列里面的格式决定。
输出文件名为tran_8.txt - bedtools转换参数
bedtools getfasta -fi /disks/backup/chaim/maize/Zea_mays.B73_RefGen_v4.42.fa -bed tran_8.txt -split -name -fo trans_fa
最后输出文件为trans_fa
4.提取我们需要的转录本--合并后的转录本
sed -n '/trans/{N;p}' trans_fa >transcript
transcript即为我们的最终结果。
