私のgtf
ファイルは100以上のディレクトリにあります。以下は、彼らがどのように見えるかを示しています。
SampleA
|___________ SampleA.GRCh38.gtf
SampleB
|___________ SampleB.GRCh38.gtf
gtf
ここでは、2つのファイルのみを例として示します。
SampleA.GRCh38.gtf
次のように:
# stringtie -e -B -p 8 -G /path/stringtie_output/stringtie_merged.gtf -o /path/SampleA.GRCh38.gtf /path/SampleA.sorted.bam
# StringTie version 1.3.3
chr1 StringTie transcript 11594 191502 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.2"; cov "0.0"; FPKM "0.000000"; TPM "0.000000";
chr1 StringTie exon 11594 14829 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.2"; exon_number "1"; cov "0.0";
chr1 StringTie exon 14970 15038 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.2"; exon_number "2"; cov "0.0";
chr1 StringTie exon 15796 16765 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.2"; exon_number "3"; cov "0.0";
chr1 StringTie exon 16858 17055 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.2"; exon_number "4"; cov "0.0";
chr1 StringTie exon 17233 17742 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.2"; exon_number "5"; cov "0.0";
chr1 StringTie exon 17915 18061 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.2"; exon_number "6"; cov "0.0";
chr1 StringTie exon 18268 19364 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.2"; exon_number "7"; cov "0.0";
chr1 StringTie exon 189836 191502 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.2"; exon_number "8"; cov "0.0";
chr1 StringTie transcript 11594 195411 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.6"; cov "0.0"; FPKM "0.000000"; TPM "0.000000";
chr1 StringTie exon 11594 14829 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.6"; exon_number "1"; cov "0.0";
chr1 StringTie exon 14970 15236 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.6"; exon_number "2"; cov "0.0";
chr1 StringTie exon 185758 187287 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.6"; exon_number "3"; cov "0.0";
chr1 StringTie exon 187376 187577 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.6"; exon_number "4"; cov "0.0";
chr1 StringTie exon 187755 187890 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.6"; exon_number "5"; cov "0.0";
chr1 StringTie exon 188130 188266 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.6"; exon_number "6"; cov "0.0";
chr1 StringTie exon 188439 188584 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.6"; exon_number "7"; cov "0.0";
chr1 StringTie exon 188791 188902 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.6"; exon_number "8"; cov "0.0";
chr1 StringTie exon 195263 195411 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.6"; exon_number "9"; cov "0.0";
chr1 StringTie transcript 11594 197912 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.5"; cov "0.0"; FPKM "0.000000"; TPM "0.000000";
次のSampleB.GRCh38.gtf
ように:
# stringtie -e -B -p 8 -G /path/stringtie_output/stringtie_merged.gtf -o /path/SampleB.GRCh38.gtf /path/SampleB.sorted.bam
# StringTie version 1.3.3
chr1 StringTie transcript 11594 191502 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.2"; cov "0.0"; FPKM "0.000000"; TPM "1.000000";
chr1 StringTie exon 11594 14829 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.2"; exon_number "1"; cov "0.0";
chr1 StringTie exon 14970 15038 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.2"; exon_number "2"; cov "0.0";
chr1 StringTie exon 15796 16765 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.2"; exon_number "3"; cov "0.0";
chr1 StringTie exon 16858 17055 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.2"; exon_number "4"; cov "0.0";
chr1 StringTie exon 17233 17742 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.2"; exon_number "5"; cov "0.0";
chr1 StringTie exon 17915 18061 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.2"; exon_number "6"; cov "0.0";
chr1 StringTie exon 18268 19364 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.2"; exon_number "7"; cov "0.0";
chr1 StringTie exon 189836 191502 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.2"; exon_number "8"; cov "0.0";
chr1 StringTie transcript 11594 195411 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.6"; cov "0.0"; FPKM "0.000000"; TPM "3.000000";
chr1 StringTie exon 11594 14829 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.6"; exon_number "1"; cov "0.0";
chr1 StringTie exon 14970 15236 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.6"; exon_number "2"; cov "0.0";
chr1 StringTie exon 185758 187287 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.6"; exon_number "3"; cov "0.0";
chr1 StringTie exon 187376 187577 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.6"; exon_number "4"; cov "0.0";
chr1 StringTie exon 187755 187890 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.6"; exon_number "5"; cov "0.0";
chr1 StringTie exon 188130 188266 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.6"; exon_number "6"; cov "0.0";
chr1 StringTie exon 188439 188584 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.6"; exon_number "7"; cov "0.0";
chr1 StringTie exon 188791 188902 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.6"; exon_number "8"; cov "0.0";
chr1 StringTie exon 195263 195411 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.6"; exon_number "9"; cov "0.0";
chr1 StringTie transcript 11594 197912 . - . gene_id "MSTRG.7542"; transcript_id "MSTRG.7542.5"; cov "0.0"; FPKM "0.000000"; TPM "0.000000";
transcript
3番目の列からのみ抽出したいです。transcript_id
どの列が10番目の列で、どの列がTPM
最後の列かを抽出したいと思います。ただし、TPM
サンプル名でなければなりません。
出力が次のようになります。
Type transcript_id SampleA SampleB
transcript MSTRG.7542.2 0.000000 1.000000
transcript MSTRG.7542.6 0.000000 3.000000
transcript MSTRG.7542.5 0.000000 1.000000
答え1
各ファイルから関連レコードを抽出し、結果を2つの新しい一時ファイルに書き込む必要があります(おそらくawk
キーソートsort
を使用して)。以下は、ファイルの1つを処理する例です。
awk '$3 == "transcript" {printf("%s %s %s ", $3, $10, $12, $18);}' SampleA.GRCh38.gtf | sort -k 2 > tf1
その後、各レコードに各ファイルの最後の2つの列が含まれるようにjoin
マージして作成された2つの一時/中間ファイルを使用できます。awk
join
使用できるコマンドの例は次のとおりです。
join -o 1.1,1.2,1.3,2.3 -1 2 -2 2 tf1 tf2
実行する前に、ヘッダー行を印刷して(例:コマンドを使用printf
)、出力のスペースをタブに置き換える(または使用する)、他のスクリプトを使用して出力フォーマットを指定join
できます。join
sed
awk
これらの例では、両方のファイルを処理し、目的の出力を生成する(そして一時ファイルをクリーンアップするなど)スクリプトを作成できるはずです。
データファイルのサイズに応じて、1つまたは複数のプログラムでawk
すべての操作を実行できます(つまり、両方のファイルで選択したすべてのデータを同時にメモリに簡単に保存できます)。python
perl
答え2
関心のある行にのみ18番目のフィールドがあるので、ファイルを削除してからjoin
そのawk
ファイルを含むファイルを削除するだけです。NF==4
他のすべての行には2つのフィールドしかありません。
計算されたパスの特定の仮定もありますが、SampleB
適切に変更できます。
while IFS= read -r -d '' f; do #read the list of SampleA
g=$(echo "$f" | sed "s/pleA/pleB/g") #calculate path to SampleB
if [[ -f "$g" ]]; then #check SampleB exists
echo "$f" | sed "s/.*pleA\.//g" #print sample No
echo "Type transcript_id SampleA SampleB" #print header
#do the join
join -j 12 -o 1.3 -o 1.12 -o 1.18 -o2.18 <(sort -k 12 "$f") <(sort -k 12 "$g") | awk 'NF==4'
fi | sed 's/[;"]//g'| column -t #make it pretty
done < <(find . -type f -iname "*SampleA*" -print0) #NULL separated list of SampleA
答え3
次のコマンドを試してください
ステップ1
awk '$3 ~ /transcript/{print $0}' file1|awk '{print $3,substr($12,2,12),substr($NF,2,8)}' > out1
ステップ2
awk '$3 == "transcript" {print substr($NF,2,8)}' file2 > out2
ステップ3
paste out out1.txt | awk 'BEGIN{print "Type transcript_id SampleA SampleB"}{print $0}'
Output
Type transcript_id SampleA SampleB
transcript MSTRG.7542.2 0.000000 1.000000
transcript MSTRG.7542.6 0.000000 3.000000
transcript MSTRG.7542.5 0.000000 0.000000